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ABSTRACT 


The  second-order  solution  of  the  problem  of  the  inter- 
action of  a train  of  regular  waves  with  a completely  sub- 
merged, horizontal  circular  cylinder  in  finite  depth  water 
is  presented  for  the  two-dimensional  case.  The  incident 
wave  is  developed  as  a second-order  Stokes'  wave  by  use  of 
a perturbation  method  and  the  solution  of  both  the  first- 
order  and  second-order  scattering  potentials  is  obtained 
numerically  using  the  Green's  function  approach.  The  hydro- 
dynamic  pressure  and  resulting  first-order  and  second-order 
force  coefficients  are  determined  numerically  and  presented 


for  various  values  of  water  depth,  cylinder  depth  of 
submergence,  and  wave  length. 
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Description 
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Subscripts 
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Symbol 


Description 
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The  solution  to  the  two-dimensional  linear  wave/structure 
interaction  problem  has  been  well-studied  for  both  the  finite 
and  infinite  depth  cases.  The  fluid  motion  resulting  from 
the  interaction  of  a train  of  regular  waves  with  a submerged 
horizontal  circular  cylinder  in  infinite  depth  water  was 
first  studied  by  Dean  [Ref.  1] . Ursell  [Ref.  7]  later  studied 
the  problem  anew,  placing  the  solution  on  a rigorous  basis. 

More  recently,  Ogilvie  [Ref.  6]  reconsidered  the  problem.  He 
computed  the  first-order  oscillatory  forces  and  second-order 
steady-state  forces  for  the  following  cases:  (a)  the  cylinder 

held  fixed,  (b)  the  cylinder  forced  to  oscillate  sinusoidally 
in  still  water,  and  (c)  the  neutrally  buoyant  cylinder  allowed 
to  respond  to  wave  motion. 

It  can  easily  be  shown  that  the  steady-state  part  of  the 
second-order  forces  can  be  computed  from  the  first-order 
potential  alone.  Accordingly,  Ogilvie  did  not  obtain  a 
complete  second-order  solution;  the  steady-state  forces 
arise  from  the  first-order  velocity- squared  terms  in  Ber- 
noulli's equation.  In  addition  to  the  steady-state  forces, 
second-order  oscillatory  forces  also  exist  which  have  a 
frequency  twice  that  of  the  first-order  fox  *es  and  represent 
the  subject  of  the  present  work. 

This  report  reconsiders  again  the  submerged  horizontal 
cylinder  problem,  extending  the  analysis  to  include  water 

I 
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of  finite  depth.  The  cylinder  is  considered  to  be  completely 
submerged  and  held  fixed  in  a train  of  regular  waves.  The 
objective  is  to  determine  the  complete  second-order  solution 
and,  accordingly,  determine  the  oscillatory  second-order 
forces  as  well  as  the  steady-state  second-order  forces. 

The  problem  is  treated  as  a regular  perturbation  problem 
in  the  small  parameter,  incident  wave  height/ cylinder  radius. 
Proceeding  in  this  way  the  incident  wave  appears  as  a second- 
order  Stoke' s wave.  The  boundary-value  problems  for  both 
the  first-order  and  second-order  potentials  are  established 
and  the  solutions  to  both  are  obtained  by  use  of  the  Green's 
function  method. 
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II.  THEORETICAL  DEVELOPMENT 


A.  FORMULATION  OF  THE  PROBLEM 

A systematic  representation  of  the  problem  considered  is 
illustrated  in  Fig.  (1) . A rigid  right  cylinder  of  radius  a, 
submerged  to  a depth  d in  water  of  depth  h,  is  fixed  in 
a train  of  regular  waves  propagating  in  the  positive  x direc- 
tion. The  basic  problem  is  that  of  calculating  the  hydro- 
dynamic  pressure  on  the  immersed  cylinder  and  the  resulting 
force  correct  to  the  second-order  in  the  wave  height  of  the 
incident  wave. 

Assuming  the  fluid  to  be  irrotational , a velocity 
potential,  <5,  may  be  defined  as: 

q = V$(x,y,t)  (1) 

where  q denotes  the  fluid  velocity  vector.  (The  barred 
quantities  denote  dimensional  quantities.)  Moreover,  assuming 
the  fluid  to  be  incompressible,  and  homogeneous,  the  velocity 
potential  must  satisfy  the  Laplace  equation, 

V2$(x,y,t)  = 0 (2) 

within  the  fluid  region. 

In  addition  to  the  differential  equation,  Eq.  (2) , $ 
must  satisfy  certain  boundary  conditions.  In  specific, 
these  are: 
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1.  The  velocity  of  the  fluid  normal  to  the  bottom 


boundary  must  vanish. 


2.  The  velocity  normal  to  the  surface  of  the  cylinder, 


defined  by  S^(x,y)=0,  must  vanish. 


3.  The  kinematic  boundary  condition  on  the  free  surface. 


defined  by  S2  (x,y , t) = o (x, t) =0 , must  be  satisfied. 


4.  The  dynamic  boundary  condition  on  the  free  surface 


must  be  satisfied. 


These  boundary  conditions  may  be  stated  mathematically  as  follows: 


4>-(x,-h,t)=0 


V4>-VS1=0 


q-S-  + 


2 3 1 


$-(x,n,t)+  *5  [<t>- (x,  n,  t)  + $-(x,n,t)  ]+gn(x,t)=gK  (6) 


where  n denotes  the  elevation  of  the  free  surface,  g denotes  the 


acceleration  of  gravity,  and  K denotes  the  Bernoulli  constant. 


For  convenience  in  carrying  out  the  solution,  the  variables 


are  next  made  dimensionless  using  the  cylinder  radius,  and  the 


wave  frequency,  a , as  follows: 


I 


p 


1 


x/a 

y = y/a 

d = d/a 

h = h/a 

H/2a 

K = K/a 

2— , 

v = a a/g 

t = at 

aO/ga 

n = n/a 

(7) 


H denotes  the  dimensionless  wave  height,  K denotes  the 
dimensionless  Bernoulli  constant,  and  t denotes  the 
dimensionless  time. 

Utilizing  the  parameters  defined  in  Eq.  (7),  the  boundary 
value  problem  defined  in  Eqs.  (2-6)  may  be  rewritten  concisely 
in  dimensionless  form  as: 


V <|>(x,y,t)  = 0 


in  the  fluid  region 


(8) 


$y(x,-h,t)  = 0 (9) 

<t>n(x,y,t)  =0  on  S^x^)  = 0 (10) 

$x(x,n ,t) nx(x,t)  - $ (x,n»t)  + vnt(x,t)  = o (11) 

<j>t  (x,  n » t)  + ^Ux(x,n,t) 2 + 4>y  (x,n,t)  2]  + n(x,t)  = k (12) 


B.  PERTURBATION  EXPANSION 

Expanding  <J>  and  n in  terms  of  a perturbation  parameter, 
e,  provides  a means  for  converting  the  nonlinear  boundary-value 
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problem  into  a series  of  linear  problems.  The  solution  to 
each  linear  problem  may  be  tractable  whereas  the  nonlinear 
problem  is  unsolvable. 

Since  <+> , n,  and  K are  functions  of  the  small  parameter, 
e,  they  may  be  written  in  the  power  series: 


<J>(x,y,t;e)  = £ e%n  (x,y,  t)  (13) 

n=l 


r)  (x,  t;  e)  = £ enn  (x,t)  (14) 

n=l  n 


and 


R=EenKn  (15) 


In  Eqs.  (11)  and  (12),  4>  contains  ru  and,  therefore,  c, 
implicitly.  This  can  be  converted  to  an  explicit  form  in 
e by  use  of  the  Taylor  series  expansion: 


4>(x,n,t) 


00  , . \ m 

V h (x,t ; e) 

m=0  ml 


r3md>  (x,y  ,t)  . 

m Jy=0 

3y  * 


(16) 


The  perturbation  parameter,  e,  is  related  to  the  wave  height 
in  an,  as  yet,  unknown  manner. 

Equations  (13-16)  as  well  as  appropriate  derivatives 
similar  to  Eqs.  (13-16)  may  now  be  substituted  into  the 
boundary-value  problem  given  in  Eqs.  (8-12)  to  obtain  a 
series  of  linear  boundary-value  problems  for  $2,  etc. 
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That  is,  upon  substitution  and  equating  coefficients  of  like 
powers  of  e,  the  following  problems  for  the  first  two  terms 
in  the  expansion  for  $ can  be  precipitated: 

First-order  boundary-value  problem  (e): 


V 2 (x,y,t)  = 0 (17) 

<j>ly(x,-h,t)  = 0 (18) 

♦ln(x,y,t>  = 0 on  Sx(x,y)  = 0 (19) 

4>ly(x,0,t)  - vnlt(x,t)  = 0 (20) 

n^Xjt)  + (Jlt(x,0,t)  = 0 (21) 


2 

Second-order  boundary-value  problem  (e  ): 

V2<}>2(x,y,t)  = 0 (22) 

<J>2y(x,-h,t)  = 0 (23) 

4>2n(x,y,t)  =0  on  S1(x,y)  = 0 (24) 

$2y(x,0,t)  ■ vn2t(x,t)  = 

(25) 

-n1  (x,t)  <^lyy (x,0 , t)  + 4>lx(x,0,t)nlx(x,t) 


ri2(x,t)  + 4>2t(x,0,t)  = -n1  (x,t)  (j>lyt  (x,0,t)  - 


There  is  a striking  similarity  between  the  first-order 
and  second-order  problems;  only  the  right  hand  side  of  the 
free  surface  boundary  conditions  differ.  In  the  first-order 
problem  the  free  surface  boundary  conditions  are  homogeneous 
whereas  the  second-order  free  surface  boundary  conditions 
are  non-homogeneous , the  right  hand  side  being  functions  of 
the  first-order  velocity  potential,  first-order  free  surface 
elevation,  and  their  derivatives. 

The  first-order  potential  may  be  represented  by  a function 
which  is  periodic  with  the  fundamental  frequency.  Accordingly, 
the  complex  potential,  u1(x,y),  may  be  expressed  as: 

4>i  (x ,y , t)  = abRe  [iu1  (x ,y)  e~it]  (27) 

where  Rg  denote  the  real  part,  a = 2ira/L,L  denoting  the  wave 
length,  and  b is  an  unknown  real  constant. 

C.  THE  FIRST-ORDER  BOUNDARY  VALUE  PROBLEM 

Since  the  first-order  boundary-value  problem  is  linear,  the 
total  first-order  potential  may  be  expressed  as  the  sum, 


I S 

where  denotes  the  incident  wa”e  potential  and  denotes 

the  scatter  potential  which  is  due  to  the  presence  of  the  cylinder. 

The  space  dependent  part  of  the  complex  potentials  are  expressed, 

in  view  of  Eqs . (27)  and  (28)  as. 
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ux  = + uf  (29) 

The  potential,  u*,  represents  only  a regular  wave  propa- 
gating along  the  x-axis  and  this  potential  must  satisfy  the 
first-order  boundary-value  problem  when  no  rigid  body  is 
present.  Therefore,  satisfying  Eqs.  (17-18)  and  (20-21)  , 
is  given  by 

I _ 1 cosh  [a(y+h)l  iax 

U1  a cosh  (ah)  e 

The  parameters  a and  v are  related  through  the  well-known  equation 
from  linear  wave  theory. 

2- 

v = ■■  = a tanh(ah)  (31) 

which,  in  dimensional  form,  using  a = ka , and  k = 2*/L  is: 

a^=  gk  tank(kh)  (32) 

Since  the  solution  for  u*  is  known,  the  boundary-value 

s 

problem  for  u^  may  now  be  established.  Substituting  Eqs.  (30), 

(29)  and  (27)  into  the  first-order  problem  given  by  Eqs.  (17-21), 
and  eliminating  between  Eqs.  (20)  and  (21)  results  in  a 


(30) 
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boundary-value  problem  for  the  first-order  scatter  potential 
as  follows: 


(x,y) 


V2u^(x,y)  = 0 
u^y (x, -h)  - 0 

cosh1|ah)  [ny  sinhta(y+h))  + inx 
on  S1 (x,y)  =0 

u^y (x, 0)  - vu^(x,0)  = 0 


(33) 


(34) 


cosh [a (y+h) 


lax 


(35) 


(36) 


where  nx  and  ny  denote  the  x-  and  y-components , respectively, 

A A A 

of  the  unit  normal  vector,  n = in^  + 3ny'  directed  outward 
into  the  fluid. 


D.  SECOND-ORDER  BOUNDARY-VALUE  PROBLEM 

In  view  of  the  linearity  of  the  second-order  problem, 
the  same  procedure  as  used  in  the  first-order  problem  may 
be  applied  and  leads  to  the  representation  of  the  second-order 
potential  as  the  sum: 


<t>2  = 4>2  + $2  (37) 

where  ^ denotes  the  second-order  incident  wave  potential 

S 

and  $2  denotes  the  scatter  potential. 
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Considering  the  case  where  there  is  no  body  present, 

s s 

there  is  no  scattered  wave  and  therefore,  <J>^  = <f>2  = 0 . 
Additionally,  the  boundary  conditions  on  the  immersed  cylinder, 
Eqs . (19)  and  (24),  are  not  applicable.  Thus,  when  the  substitu- 
tions of  <+>*  for  <f)^  and  f°r  4>2  are  ma^e  into  Eqs.  (22)  , 

(23),  (25)  and  (26),  along  with  Eqs.  (20),  (21),  (27)  and 
(31),  and  upon  eliminating  and  between  Eqs.  (25)  and 
(26) , the  resulting  boundary-value  problem  for  the  second- 
order  incident  wave  potential  becomes: 

V2 (x,y , t)  = 0 (38) 

(x,-h,t)  = 0 (39) 

<J>2y(x,0,t)  + V(f>2tt  (x,0  ,t)  = - jb2  (a2-v2)  Re  [iei2  (ax_t)  ] 

(40) 

The  periodic  solution  to  the  incident-wave  boundary-value 
problem  specified  by  Eqs.  (38-40)  is  known  and  may  be  expressed 
as 


$2  = |ab2vRe[iu2(x,y)e"l2t] 


(41) 


where  the  second-order  complex  potential,  U2 , is  given  by: 


u2(x,y)  = 


1 cosh [2a (y+h) ] „i2ax 
__  : _ e 

sinh  (ah) 


(42) 
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Eqs.  (41-42)  axe  familiar  in  wave  theory  and  are  exactly  the 
same  form  as  that  given  by  Ref.  4 for  second-order  Stokes 
waves . 

Having  developed  a solution  for  the  second-order  incident 
wave  potential,  the  problem  for  second-order  scatter 

I s 

potential  may  be  determined.  To  this  end,<j>^=$i  + ^ and 
4, 2= 4> 2 + <t>2  sre  substituted  into  Eqs.  (22-26),  and  the  pre- 

viously determined  solutions  for  the  incident  wave  potentials, 

Eos.  (30)  and  (42),  are  utilized.  After  applying  Eqs.  (20),  (21), 
(27),  (29)  and  (41)  and  eliminating  n2  between  Eqs.  (25)  and  (26), 

the  resulting  boundary-value  problem  for  the  second-order  scatter 
potential  becomes: 


(x,y,t) 


V2<J>2  (x,y,t)  = 0 


(43) 


<t>2y  (x,-h,t)  = 0 


(44) 


3 

T 


ab  R_  r 


sinh  (ah) 


(n  cosh [2a (y+h) ] 


- in^  sinh [2a (y+h) ])  e 


i2 (ax-t) 


(45) 


on  S1(x,y)  = 0 


¥ - ««  ««> 

1 S I .IS  - S 2 , S 2.  -i2t. 

+ — u,  u.  - 4u.  u.  - 2un  - 3u.  )e  j 

v ly  lyy  lx  lx  lx  ly 


+ U(x) 


on  y = 0 
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where  U(x)  is  a non-periodic  function  generated  by  the 


substitution  of  into  Eqs . (25)  and  (26).  For  brevity, 

this  function  is  not  written  out  since  it  is  not  needed  in 
determining  the  second-order  pressures  and  forces. 

S 

The  boundary-value  problem  in  <t>2  given  in  Eqs,  (43-46) 
is  time  dependent  according  to  e l2t,  except  for  the  U(x) 

S 

term  in  Eq.  (46).  Therefore,  the  solution  for  $ . may  be 
taken  in  the  form 


<J>2  = j ab2 vRg £i  [u^  (x,y)e  l2t  + u^^y)]]  (47) 

where  the  last  term  denotes  the  time-independent  portion  of 

the  complex  potential.  Separate  boundary-value  problems  for 
S ~ S 

U2  and  U2  arise  from  the  substitution  of  Eq.  (47)  into  Eqs. 

S 

(43-46).  However,  as  will  be  demonstrated  only  the  $2t  term 
is  required  for  determination  of  the  pressure  to  the  second 
order.  Therefore,  the  time-derivative  of  the  time-independent 

S 

part  of  <p 2 will  be  zero,  and,  accordingly,  will  not  be  developed 

S 

further  here;  the  solution  for  u2  only  will  be  pursued. 

When  Eq.  (47)  is  substituted  into  Eqs.  (43-46),  the  result- 
ing boundary-value  problem  for  the  second-order  scatter  potential 
takes  the  form: 

V2U2 (x , y)  = 0 (43) 


(x , -h) = 0 


(49) 
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U?  (x,y)  = fn  sinh  [2a  (y+h)  ] 

sinh  (ah)  L Y 

+ inx  cosh  [2a (y+h) ]J  el2ax  on  Sx(x,y)=0  (50) 

u2y(x»0)  - 4vu2(x,0)  = f* (x)  (51) 


where 


fMx) 


2arl  SS  IS  -I  S , 1 S I 

— [— u,  u,  + u,u,  - 6u,  u,  + -u,  u, 

3v  v ly  lyy  1 lyy  ly  ly  v ly  lyy 


- i s 
- 4u.  u. 
lx  lx 


y=0 


(52) 


There  is  considerable  similarity  in  the  first-order 
and  second-order  problems,  the  only  difference  in  form  being 
that  the  second-order  free  surface  boundary  condition, 

Eq.  (51),  is  non-homogeneous . 

Because  of  the  non-homogenity  of  Eq.  (51),  further  use  is 
made  of  linear  superposition  defining  u^  as  the  sum, 


(53) 


S S * 

where  u^  denotes  the  homogeneous  solution  and  u^  the 

particular  solution  to  the  boundary-value  problem  as  stated 

S * S° 

in  Eqs . (48-51).  More  precisely,  U2  and  U2  are  defined  as 
the  solutions  to  the  boundary-value  problems  obtained  from 
the  substitution  of  Eq.  (53)  into  Eqs.  (43-51)  as  follows: 
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Particular  Solution 


) : 


V2U2*(x,y)  = 0 


(54) 


u^y  (x,-h)  = 0 

u^* (x,  0)  - 4vu^* (x,0)  = f * (x) 

S° 

Homogeneous  Solution  (u^  ) 

2 S° 

V u2  (x,y)  = 0 
s° 

u2y (X'"h)  = 0 

qO  gO 

u2y(x,0)  - 4vu2  (x ,0)  = 0 

u^  (x,y)  = 3 In  sinh  (2a  (y+h)  ] 

2n  sinh^(ah)  LY 

* 

+ in^  cosh [2a (y+h) ] Jei2ax  - u2n^x'^^ 
on  S1(x,y)  = 0 


(55) 


(56) 


(57) 


(58) 

(59) 


(60) 


By  the  division  of  the  second-order  problem  into  homogeneous 

S* 

and  particular  solutions,  the  non-homogeneous  problem  for  u2 

contains  no  boundary  condition  on  the  cylinder  surface.  The  result- 

S* 

ing  boundary-value  problem  for  u2  , as  stated  in  Eqs . (54-56),  is 

identical  to  that  associated  with  the  linear  problem  for  the  potential 
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resulting  from  a harmonic  pressure  variation  of  amplitude 


distribution  f*(x)  on  the  free  surface  in  water  of  depth  h. 

S° 

The  homogeneous  boundary-value  problem  for  u^  , defined  by 

Eqs . (57-60),  is  similar  in  form  to  the  first-order  problem 

given  by  Eqs.  (33-36),  the  significant  difference  being  the 

term  4v  in  place  of  v which  occurs  in  the  free  surface  boun- 

s° 

dary  condition.  Thus,  the  method  of  solution  for  u2  will 
be  similar  to  that  of  the  first-order  scatter  potential 

E.  DEFINITION  OF  THE  INCIDENT  WAVE 

Having  developed  the  first-order  and  second-order  boun- 
dary-value problems,  it  is  now  appropriate  to  completely 
define  the  incident  wave  height  and  in  the  process  determine 
expressions  for  the  unknown  constants  b and  . Solving 
Eqs.  (21)  and  (26)  for  and  respectively,  and  then 
substituting  the  results  into  the  expression  for  the  free 
surface  elevation  as  given  by  Eq.  (14)  yields: 


n(x,t)  = - c*lt  + e [-*2t  + *lt*lyt  + ♦ltt*ly 


(61) 


- + *!y>  + K2]  + 0 ( e3) 


where  the  potentials  are  evaluated  at  y = 0 . Eq.  (61) 
expresses  the  free  surface  elevation  in  terms  of  the  total 
potential,  and  therefore,  includes  the  effects  of  both  the 
incident  and  the  scattered  waves.  Hence,  as  it  is  of  inter- 
est to  obtain  an  expression  for  the  free  surface  elevation 


of  the  incident  wave,  the  scatter  potentials  are  stt  to 

s s 

zero,  i.e.  ^ = 0 . Evaluating  Eq.  (61)  using  the  known 

solutions  for  the  incident  wave  potentials,  Eqs . (30)  and 
(42),  along  with  Eqs.  (27)  and  (41),  then  yields: 


(x,t>  = «bR.l.1I“-«l  * e2[b-2.<-^  ♦ 


K. 


(62) 


aD  cosh  (ah)  [2  + cosh  (2ah)  ] n r_i2 (ax-t) , 

+ —j—  = e ie  J 

sinh-3  (ah) 


+ 0 (e  ) 

where  nI(x,t)  denotes  the  dimensionless  surface  elevation 
for  the  incident  wave  with  no  cylinder  present.  If  the  x-axis 
is  placed  at  the  mean  water  line  then  the  constant  term  must 
vanish  and,  therefore. 


u2,  2 2, 

„ _ b (a  - v ) 

K2 TS 


(63) 


The  remaining  terms  in  Eq.  (62)  are  time -dependent,  periodic 
functions,  the  second-order  term  having  a frequency  of  twice 
that  of  the  first-order. 

Defining  the  dimensionless  wave  height  given  in  Eq.  (7) 
as  the  difference  in  elevation  between  the  crest  and  trough 
of  the  incident  wave,  eb  must  represent  the  wave  amplitude. 
The  second-order  term  in  Eq.  (62)  having  frequency  twice  the 
fundamental  frequency  makes  equal  contributions  to  surface 
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elevation  at  both  the  crest  and  trough,  so  that  it  contributes 
nothing  to  the  wave  height.  Thus,  in  terms  of  dimensionless 
wave  height,  b is  given  by 


H = eb  = H/2a  (64) 

where  H denotes  the  elevation  difference  between  the  wave 
crest  and  trough. 

Expressing  the  incident  wave  profile  in  terms  of  the 
dimensionless  wave  height  yields: 

nI(x,t)  = H cos  (ax-t)  (65) 

+ H.  a-E£s^-(?h.)  [2  + cosh  (ah)  ] cos  [2  (ax-t)] 

4 sinhJ(ah) 

It  may  be  noted  that  Eq.  (65)  agrees  with  the  expression 
for  the  second-order  Stokes'  wave  given,  for  example,  by 
Ref.  4. 

F.  PRESSURES  AND  FORCES 

The  pressure  may  be  formulated  by  use  of  Bernoulli's 
equation  arranged  as  follows: 

2 2 

P(x,y,t)  = - jp[$-  + ] - pgy  + pgK  (66) 

where  P(x,y,t)  denotes  the  dimensional  pressure  and  p denotes 
the  fluid  density.  By  use  of  Eqs . (7)  , (13-16)  , (27)  , (41)  , 
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(47) , and  (64) , Eq.  (66)  may  be  reduced  to  dimensionless  form, 
and  carried  to  the  second-order  in  wave  height.  The  result  is 


P(x,y,t)  = -y  - HaRe[u1e”lt] 


(67) 


,2  2 


6v“  2 2,  -i2t, 

e u2  - ulx  ' uly>  ‘ 


H ^ n rUV  \ j-  V,  «• 

- - u,„)e  ] 


+ lulx(  + |uly 


2 , 4 - y] 

a J 


In  Eq.  (67)  p(x,y,t)  denotes  the  dimensionless  pressure 
coefficient  and  is  defined  as: 


p (x,y , t)  = P-(Xf^-  (68) 

pga 

The  first  term  in  Eq.  (67)  represents  the  hydrostatic 
pressure  as  y is  the  dimensionless  depth  beneath  the  mean 
free  surface  (y  = 0) . The  second  and  third  terms  are  harmonic, 
the  third  term  having  twice  the  frequency  of  the  second,  and 
representing  the  harmonic  second-order  contribution.  The 
remaining  terms  in  Eq.  (67)  are  independent  of  time  and  result 
in  the  time-average  or  steady-state  force  components. 

Expressing  the  components  of  the  wave  force  vectors  in 
terms  of  integrals  of  the  pressure  over  the  cylinder  surface 
area,  the  dimensionless  force  coefficients  may  be  written 
as 


ci(t) 


- Jp  (x,y , t) ni  dS1,  i = 1,2 
S1 


(69) 
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where  the  dimensionless  force  coefficients  in  the  x and  y 
directions  are  defined,  respectively,  as 


cx  (t) 


Fx(t) 

pga3 


C2  (t) 


yt) 

pga3 


(70) 


(71) 


where  F^(t)  and  F^(t)  denote  the  force  components.  Addi- 
tionally, dS-^  denotes  a dimensionless  differential  arc 
length  along  the  circumference  of  the  cylinder. 

Applying  Eq.  (67)  to  Eq.  (69)  yields: 


C±(t)  = - y fnidSi  ~ Ha  y*Re[uie  lt]ni  dsi 


(72] 


,2,3' 


. 6v 


- - uix  - v 

bi 


2 . _-i2t 


2 2 2 

+ lul*l  + lulyl  + ~ HnidS1  + 0 (H3) 

* a 


The  first  term  in  Eqi72)  results  from  the  hydrostatic 
pressure  on  the  cylinder  so  it  represents  simply  the  buoyancy 
force.  Disregarding  the  hydrostatic  pressure  the  hydrodynamic 
force  may  be  written  as 
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(73) 


Ci(t)  = H Fxi  cos(6li-t)  + H2 

+ H2[F2i  cos  ( S2i-2t)  + F“]  + 0 (H3) 


where  the  first-order,  second-order,  and  steady-state  force 
coefficients  and  phase  shift  angles  are  defined  by  comparison 
of  Eqs . (72)  and  (73)  as  follows: 


ifi 


Flie 


li 


= a 


/uinids] 


iS 


F 


2i 


2i 


2 

FSS  = 

2i  4v 


) n . dS 
l 


1 


1) nidS1 


(74) 


(75) 


(76) 


The  dimensionless  force  coefficients,  F^  and  F^^  are  real. 

The  first  term  in  Eq.  (73)  represents  the  first-order 
forces  of  the  fundamental  frequency,  a.  The  periodic  portion 
of  the  second  term  in  Eq.  (73)  represents  the  second-order 
contribution  to  the  force  at  twice  the  fundamental  frequency. 
The  last  second-order  term  represents  the  steady-state  or 
time-independent  contribution,  sometimes  called  drift-force. 

Evaluation  of  the  force  coefficients  defined  by  Eqs. 
(74-76)  is  the  main  objective  of  this  work.  Therefore, 
it  is  necessary  to  evaluate  the  first-order  and  second-order 
complex  potentials,  u^  and  u as  well  as  derivatives  of  u^ 
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on  the  surface  of  the  cylinder.  Additionally,  evaluation 
of  u^  and  its  derivatives  on  the  mean  free  surface  (y  - 0) 
is  required  in  order  to  carry  out  the  solutions. 
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III.  METHOD  OF  SOLUTION  AND  NUMERICAL  PROCEDURES 


A.  SOLUTION  OF  THE  FIRST-ORDER  PROBLEM 

The  use  of  a Green's  function  to  express  the  solution 
to  the  first-order  boundary-value  problem  was  first  formu- 
lated by  John  [Ref.  5],  and  applied  to  submerged  ellipsoids 
by  Garrison  and  Rao  [Ref.  3].  This  method  is  considered  to 
be  applicable,  in  principle,  to  arbitrary  shapes  and  is 
mathematically  the  most  straightforward.  The  Green's  func- 
tion, G,  of  unit  strength  which  satisfies  Eq.  (33)  as  well 
as  the  boundary  conditions,  Eqs . (34)  and  (36)  is  given  by: 


G (x,y ; E,  ,n;  v) 


In  r 


in  r-'  + opv  f r cosh  [ y (y+h)  ]cosh[y  (n+h)  ] 

0J  L (coshyh) ( vcoshyh-ysinh  uhT 


- e 


-yh  sinh(un) sinh(uy) 
sinfi  ( yh)  . 


cos  I x—  E,  | dy 


(77) 


. 4tt 

1 2a^h+sinh (2a^h) 


cosh [a1 (y+h) ] cosh [a^ (n+h) ] cos [a^ | x-£ ] ] 


where : 


1 

[ (x-£)  2 + 

(y-n)2]2 

(78) 

l 

[ (x-£)  2 + 

(y+n)2]2 

(79) 

The  symbol,  a^,  is  defined  in  terms  of  h and  v as  the 
solution  to  the  equation 
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F (a1 ,h, v) 


0 


(80) 


where  F is  defined  by 


F(alfh,v)  = tanh(a^h)  - v 


(81) 


Comparison  of  Eqs . (31),  (80)  and  (81)  demonstrates  that  a1 

is  clearly  the  equivalent  of  a.  However,  this  will  not  be 
the  case  for  the  second-order  problem,  thus,  requiring  the  use 
of  separate  notation.  In  Eq . (77),  PV  denotes  the  principal 

value  of  the  integral. 

The  Green's  function  given  in  Eq . (77)  was  also  given  in 

series  form  by  John  [Ref.  5]  as 

2 2 i 

« ( Ui+v  )cos[y.  (y+h)  ] -U.  |X-5 

G(x,y;C,n;v)  = 2tt  E 5 3 cos[y  (n+h)]e  ' 

k=l  vyk_hVik(^+,J  ) 


iailx-5| 

4:Tcosh  [a1  (y+h)  ] cosh  [a^  (n+h)  ]e 


2a^h  + sinh (2a^h) 


(82) 


where  a.,  is  as  defined  in  Eq.  (80)  and  the  quantities  p. 
are  defined  as  the  real  positive  roots  of  the  equation 


K(uk,h,v)  = 0 

where  the  function  K is  defined  as 


(83) 
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s 

Following  the  Green's  function  method  of  solution,  u^ 
is  written  as  the  integral  over  the  cylinder  arc  length, 

S as: 


u^(x,y) 


C,n)G(x,y;C/n;v) 


dS. 


(85) 


where  (£,n)  denotes  points  on  the  immersed  surface,  f^(£,n) 
denotes  the  unknown  source  strength,  and  dS^  = dS^/a  denotes 
the  differential  arc  length  on  the  cylinder  surface,  made 


dimensionless  with  the  characteristic  length,  a. 

Although  the  solution  to  the  first-order  boundary-value 
problem  as  stated  in  Eqs . (33-36)  is  given  by  Eq . (85) , the 

source  strength  function,  f^(£,n)»  must  be  determined  in 
order  to  evaluate  the  potential.  From  potential  theory, 

S 

the  normal  derivative  of  the  potential,  u^(x,y),  on  the 


surface  of  the  cylinder  is 


um(x,y) 


= y vx'Y) 


2v  s/f  1 


(C ,n)  Gn (x,y ; C,n ;v)  dSj^  (86) 


where  Gn  (the  normal  derivative  of  G)  may  be  determined  by 
dif f ei-entiation  of  either  Eq.  (77)  or  (82)  in  a straight- 
forward manner. 

Applying  the  final  boundary  condition,  Eq . (35),  leads 

to  an  integral  equation  which  may  be  solved  for  f^, 
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2 fx(x,y)  + 2tt 


/f1(C,Tl)Gn(x,y;C,n,v)  dSL 


(37) 


cosh  (ah)  [nysinh[a(y+h)]  + in^osh  [a  (y+h)  ] 


lax 


The  solution  for  from  Eq . (87)  may  then  be  used  in  Eq. 

g 

(85)  to  determine  the  potential,  u^ . 

B.  SOLUTION  OF  THE  SECOND-ORDER  PROBLEM 

The  similarity  in  form  of  the  boundary-value  problem 
for  the  homogeneous  part  of  the  second-order  potential, 

Eqs . (57-60),  to  the  first-order  potential  is  now  utilized. 

Since  the  only  significant  difference  occurs  in  Eq.  (59) , 

S° 

where  v is  replaced  by  4v,  we  may  represent  u^  in  a form 
similar  to  Eq.  (85)  as 

u2  (x,y)  = Jf-2  (£,n)G(x,y;C,v;4v)  dS^.  (88) 

S1 

where  G (x,y; £, n ; 4v)  is  defined  by  replacing  v by  4v  in  Eqs. 
(77)  and  (82).  Therefore,  a2  is  defined  by 


F(a2,h,4v)  = 0 


(89) 


and  a^  is  replaced  by  a2  in  Eqs.  (77)  and  (82), also.  In  the 
case  of  Eq.  (82)  , Uy.  is  defined  as  the  real  positive  roots  of 
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K (uk,h,4v) 


0 


(90) 


The  source  strength  function,  f2(£,n),  appearing  in 
Eq.  (88)  is  again  determined  by  applying  the  kinematic 
boundary  condition  on  the  cylinder  surface  as  given  by  Eq . 

(59) . The  integral  equation  for  f2  may  then  be  given  by 

| f 2 (x,y)  + ± yf2(C,n)Gn(x,y;£,n;4v)  dsx  (91) 

S1 

'n  sinh  [2a  (y+h)  ] + in  cosh [2a (y+h) ] 

= y 2 

5 

sinh  (ah) 
on  S1(x,y)  = 0 

S * 

However,  to  solve  Eq.  (91)  for  f2,  first  a solution  for  u2 

S * 

must  be  obtained  and  u2n  evaluated  on  the  cylinder  surface. 

S * 

Thus,  at  this  point  the  solution  to  u2  is  sought. 

S* 

The  complex  potential  u2  is  defined  as  the  solution  to 
the  boundary-value  problem,  Eqs . (54-56).  The  form  of  this 

problem  is  recognized  as  being  equivalent  to  that  of  fluid 
motion  produced  by  a periodic  pressure  distribution  (as 
specified  by  f*(£))  with  frequency  2a  on  the  free  surface.  No 
cylinder  is  considered  to  exist  in  the  fluid  region. 

The  solution  to  this  problem  may  be  formulated  as  an 
integral  over  the  free  surface  extending  from  negative 
infinity  to  positive  infinity.  Using  the  present 


i2ax  S*  . , 

e - u2n(x,y) 
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dimensionless  representation,  U2  becomes: 

^2  (x,y)  = J?  s J** (5)G* (x,y; C,0;4v)  d £ 


The  Green's  function,  G* , for  this  problem  is  given  in  Ref.  8. 
Upon  transforming  the  solution  to  the  non-dimensional  form 


G*  (x  ,y ; 0 ; 4v)  = 


_2pv  j cosh  [ u (y+h)  ]cos[u(x-g)  Idi: 
4vcosh(ph)  - ysinh(yh) 


47Tcosh  [a2h]  cosh  [a2  (y+h)  ] cos  [a2  (x-  £)  1 


2a2h  + sinn(2a2h) 


where  a2  is  defined  by 


F(a2,h,4v)  = 0 


The  normal  derivative  of  u^  is  required  in  order  to 
completely  define  the  kinematic  boundary  condition  on  the 
cylinder  as  given  in  Eq.  (60).  From  Eq . (92),  this  becomes 


u2n(x'y)  = TF  (x,y;£,0;4v)  d£ 


The  derivative  of  G*  in  the  direction  normal  to  the  cylinder 
surface  may  be  obtained  by  differentiation  of  Eq.  (93)  in  a 
straightforward  manner.  Thus,  using  the  values  for  f*(£), 
as  defined  in  terms  of  the  first-order  solution  only, and  as 
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s * s * 

given  in  Eqs.  (52)  and  (56)  , u2  and  u2n  may  be  determined. 

s*  g o 

Using  u2n  , u2  may  then  be  determined  to  complete  the 
second-order  boundary-value  problem  solution. 

C.  NUMERICAL  METHODS 

Determination  of  the  forces  on  the  cylinder  surface 
requires  solving  for  the  potentials  u^ , u2  and  the  derivatives 
of  u-^  at  points  on  the  surface,  as  shown  in  Eq.  (72)  . Thus, 
the  problem  is  now  one  of  solving  for  both  the  first-order 
and  the  second-order  scatter  potentials,  and  requires  solving 
for  the  source  strength  functions  f^  and  f2,  as  well  as  the 
free  surface  pressure  distribution  source  strength  function, 
f*.  in  view  of  the  complexity  of  the  equations  it  is  natural 
to  attempt  a numerical  solution. 

S 

The  first-order  scattering  potential,  u^ , as  given  by 
Eq.  (85)  is  dependent  upon  the  first-order  source  strength 
function,  f ^ . Thus,  it  is  necessary  to  solve  Eq.  (87)  for 

s 

f^  to  determine  u^.  A numerical  solution  may  be  developed 
by  dividing  the  cylinder  surface  into  elements  of  length 
A 9 = 2ir/m,  with  the  center  of  each  element  assigned  an  index. 
Since  f^(£,n)  is  recognized  as  a well-behaved  function  for  a 
smooth  surface  cylinder,  it  is  appropriate  to  define  parameters 


a.  . (v)  = - 

ID  TT 


isij 


/Gn(xi'yi;C'n;v)  dSl 


i , j = 1 , 2 , 3 , . . . , m 


(96) 
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hi  ' cosh  (ah)  [yxi'yi)  sinh[a(yi+h,! 


+ m 


-1  -1.  daft.  . 

in  (x . ,y . ) cosh [a (y . +h) ] e 

X X x 


fu  ■ *x (xi'yi’ 


and, thus, Eq.  (87)  may  be  approximated  by  the  complex  matrix 
equation. 


fli  + aij(v)flj  = 2hlj 


i / j 1,2, ...,m 


Upon  evaluation  of  ci^j (v)  using  Eq . (96),  the  inversion  of 

Eq.  (99)  may  be  carried  out  on  the  digital  computer  to 
determine  f^  at  each  nodal  point  on  the  surface  of  the 
cylinder . 

For  the  purpose  of  evaluating  a and  8,  either  Eq . (77) 
or  Eq . (82)  may  be  used.  For  evaluation  of  a and  8 when  i = j 
Eq.  (77)  must  be  used  to  allow  separate  treatment  of  the 
logarithmic  singularity  as  r approaches  0.  The  difficulty 
encountered  with  the  logarithmic  singularity  may  be  overcome 
by  carrying  out  its  integration  analytically.  Garrison 
[Ref.  2]  showed  for  the  calculation  of  cu  j that  the  contribu- 
tion of  the  normal  derivative  of  the  In  (r)  in  Eq.  (77) 
integrated  over  the  singular  element  of  A0  is  A9/2,  not  only 
for  the  nodal  point  (x^y^)  = (£,n),  but  for  all  nodal  points 
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(x^,y\)  on  the  cylinder  surface.  Additionally,  to  calculate 
B — Garrison  has  developed  a numerical  approximation  for  the 
integral  of  the  In  (r)  over  the  singular  element  of  £6. 

A second  singularity  arises  in  the  evaluation  of  the 
infinite  integral  occurring  in  Eq.  (77) . The  first  term  in 
the  integrand  is  singular  at  y = a . Since  the  numerator 
approaches  zero  like  (y-a)  near  a,  Garrison  [Ref.  2]  sub- 
tracted l/(y-a)  from  the  singular  term  in  the  range  0 <_  y <_  2a 
and  added  the  contribution  of  the  integral  of  1/ (y-a) (which 
in  this  case  was  zero  ) . This  converted  the  singular  integrand 
to  a regular  function.  Applying  this  technique  to  numeri- 
cally evaluate  the  resulting  integral  between  0 and  2a, 
Simpson's  three-eights  rule  is  used  with  an  odd  number  of 
subdivisions,  thus  ensuring  that  the  point  y = a is  not 
encountered . 

With  f^  determined  from  the  inversion  of  Eq.  (99),  the 
first-order  potential  and  its  derivatives  may  be  evaluated 
on  the  cylinder  surface.  Replacing  the  surface  integrals 
with  summations,  Eq.  (85)  and  its  derivatives  may  be  written 


uli  ' 8ij  <v>  flj 


i , j 1 , 2 , 3 , . . . ,m 


(100) 


u..  . = B • • (v)  f,  . 

lxi  XI]  1] 


(101) 


ulyi  ' Byij(u,flj 


(102) 
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s s 

where  u^,  u^^,  etc.  denote  functions  evaluated  at  the 
i*"*1  nodal  point  on  the  cylindrical  surface.  The  complex 
matrices  in  Eqs . (100-102)  are  defined  as  follows: 


Bfj (v)  = ^/^(xi'yi; dSi  ifj  = 1,2, ...,m  (103) 


ASlj 


Sxij(v) 


Byij(v) 


Tn  c yGx(xi,yi;S,n;v)dS1 

ASlj 

5 = Yn  yGy  (xi,yi;C,o;v)  dS1 
ASlj 


(104) 


(105) 


The  first-order  incident  potential  and  its  derivatives  may 
be  evaluated  directly  at  each  nodal  point  after  differentia- 
tion of  Eq.  (30)  as  needed,  and  combined  with  the  scattering 
potential  and  its  respective  derivatives. 

To  solve  the  second-order  problem,  the  function  f*(£) 
required  in  Eq.  (95), and  defined  by  Eq.  (52),  must  be 
evaluated  numerically.  Dividing  the  mean  free  surface 
(y  = 0)  in  the  vicinity  of  the  cylinder  into  n equal  incre- 
ments, f*  is  evaluated  at  each  nodal  point.  The  numerical 
solution  again  uses  Eqs.  (100-105)  with  y^  = 0 for  the 

g 

purpose  of  evaluating  u^  and  its  derivatives  at  the  mean 
free  surface  nodal  points.  The  first-order  wave  potential 
and  its  derivatives  are  evaluated  directly  at  each  nodal 
point,  and  combined  with  the  scatter  potential  and  its 
derivatives  to  ^olve  Eq.  (52)  for  f*. 
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Having  determined  f*  at  all  nodal  points  on  the  free 

S* 

surface,  the  boundary-value  problem  for  u2  given  by  Eqs . 

S*  S* 

(54-56)  may  be  solved,  i.e.  u2  and  u2n  may  be  evaluated 
on  the  cylinder  by  solution  of  Eqs.  (92)  and  (95)  respec- 
tively. Expressing  the  integrals  as  numerical  summations  in 
complex  matrix  form  yields: 


S* 

u2ni 


f>a*  . ( 4v) 
3 13 


l 1 , 2 , . . . ,m 

3 1,2,.. .,n 


(106) 


ftSJj ( 4v) 


where 


(107) 


★ 

ai j ( 4 v) 
Bfj ( 4 v) 


1 

7T 


/«MW<.nMv)  dS 
2j 

57  AO  fez<xi.yiit.n,4V>  at 


(108) 


(109) 


s * s * 

Thus,  using  f * , both  u2  and  u2n  may  be  directly  evaluated 
at  each  cylinder  surface  nodal  point. 

To  solve  the  boundary-value  problem  for  the  second  part 

c ^ 

of  the  second-order  scatter  potential,  u2  , as  specified 
by  Eqs.  (57-60),  Eq.  (88)  must  be  evaluated.  However,  to 

S° 

evaluate  Eq.  (88)  for  u2  , the  second-order  source  strength 
function,  f2,  must  be  evaluated  by  numerically  solving 
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the  integral  equation,  Eq.  (91).  Replacing  the  integral 
equation  with  a complex  matrix  summation: 


f2i  + f2jaij(4v)  = 2ki  = 1'2>" -'m  (HO) 

where  a^j  (4v)  is  defined  by  Eq.  (96)  , 


1 2 j f2(xj'^j) 


and , 


k-i  = T [n  (x.  ,y . ) sinh [2a (y ,+h) ] 

L <5  inh 4 (phi  L i 1 1 1 


sinh  (ah) 


(111) 


(112) 


+ inx(xi,yi) cosh [2a (yi+h) ] 


-1  i2axi  s* 
e u2n (xi 'y,) 


Once  <*j_j(4v)  is  evaluated,  Eq.  (110)  may  be  solved  to  obtain 
the  second-order  source  strength  function,  at  each  nodal 

point  on  the  cylinder  surface.  Stating  Eq.  (88)  in  summation 
form. 


s 

u2i  = f2j8ij(4v)  i,j  = ^ » 2 , 3 , . . . ,m  (113) 

where  3.  . (4v)  is  defined  by  Eq.  (103).  Using  Eq.  (113), 

S° 

u_  may  be  determined  at  each  nodal  point  on  the  cylinder 

S* 

surface,  then  combined  with  u^  and  the  second-order  inci- 
dent potential,  u*,  to  evaluate  the  total  second-order 
potential,  u0 . 
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Since  the  first-order  potential,  u^,  and  its  derivatives, 
ulx  and  u^,  as  well  as  the  second-order  potential,  u2 , are 
now  known  at  each  nodal  point  on  the  surface,  the  first- 
order,  second-order  and  steady-state  force  coefficients  and 
phase  shift  angles  may  be  determined  using  Eqs . (74-76).  As 

shown  earlier,  the  integrals  may  be  replaced  by  summations, 
using  nodal  point  values  and  the  arc  length  increment,  £6. 

Once  each  integral  is  evaluated,  then  the  force  coefficient 
becomes  the  absolute  value  or  modulus  of  the  complex  integral 
result  and  the  phase  shift  angle  is  the  angle  whose  tangent 
is  the  imaginary  part  over  the  real  part  of  the  complex 
integral  result. 

D.  COMPUTER  SOLUTION 

In  order  to  utilize  the  method  of  solution  described 
in  Sections  B and  C,  a computer  program  was  developed  to 
carry  out  the  indicated  calculations.  Although  accuracy  was 
a prime  consideration,  an  equally  important  requirement  was 
to  minimize  computer  run  time.  To  achieve  this,  every  attempt 
was  made  to  utilize  symmetry  and  to  generate  certain  constants 
and  matrices  to  eliminate  redundant  computations. 

The  program  consists  of  a main  program  that  flows  in  the 
order  of  computation  presented  in  Sections  B and  C,  along 
with  four  subroutines.  Two  subroutines  to  solve  the  first- 
order  Green's  function,  using  both  forms  of  the  function, 

Eqs.  (77)  and  (82),  were  developed  by  Garrison  [Ref.  2]. 

These  subroutines,  GREEN  and  GREENS  respectively,  have  been 
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modified  to  include  calculation  of  the  first-order  scatter 

potential  derivatives,  evaluation  of  the  first-order  scatter 

potential  and  its  derivatives  on  the  mean  free  surface,  y = 0, 

and  evaluation  of  the  modified  Green's  function,  G*  , for 

S * S * 

the  determination  of  u^  and  U2n  at  each  cylinder  nodal 
point . 

For  elements  of  the  a and  3 matrices  corresponding  to 
small  values  of  the  parameter  (x-£) , GREEN  is  used  while  for 
larger  values  of  (x-£) , GREENS  is  used.  With  the  exception 
of  the  diagonal  elements  in  Eqs.  (96)  and  (103-105)  and  a few 
surface  nodal  points  in  Eqs.  (108-109),  the  majority  of 
elements  of  a and  8 are  calculated  by  GREENS.  This  is  most 
fortunate  as  the  series  form  converges  rapidly,  requiring  much 
less  computer  time  than  the  equivalent  integral  form. 

Subroutine  GEODAT  reads  the  input  geometrical  data, 
generates  the  matrices  h^  and  k„.  as  defined  by  Eqs.  (97)  and 
(112)  respectively,  and  calculates  certain  geometrical 
parameters  and  matrices  for  repeated  use  in  GREEN  and  GREENS. 
Subroutine  COMAT  inverts  complex  matrix  equations  and,  there- 
fore is  used  to  invert  Eqs.  (99)  and  (110)  to  determine  f^  and 
f2  respectively. 

A cross-reference  between  text  and  computer  program 
nomenclature  is  given  in  Table  I . 


TABLE  I:  Computer  Program  — Text  Symbol  Cross-Reference 
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Text 

Computer 

Program 

Text 

Computer 

Program 

a 

A 

n 

ANX(I) 

X 

d 

D 

ny 

ANY (I) 

fl 

F ( 1 , 1 ) 

ui 

U1(I) 

f 2 

Fl (1,1) 

Ulx 

U1X(I) 

f* 

FS  (L) 

uiy 

UIY  (I) 

» — 1 
rl 

u, 

Cl  (1) 

u*  (surface) 

U1IS 

F12 

Cl  ( 2 ) 

uXlx  (surface) 

UlISX 

F21 

C2  (1) 

uly  (Sur^ace) 

U1ISY 

F22 

C2  (2) 

u^yy  (surface) 

U1ISYY 

FSS 

21 

C3  ( 1) 

s 

u^  (surface) 

U1SS 

FSS 

22 

C3  ( 2) 

s 

ulx  ^sur^ace^ 

U1SSX 

G 

GI J , GI JEXT 

s 

u^  (surface) 

U1SSY 

G* 

GI J , GI JEXT 

s 

uiyy  (surface) 

UlSSYY 

G 

n 

GNIJ , GNJI 

u2 

U2  (I) 

Gx 

GXI J , GXJI 

X 

X (I) 

Gy 

GYI J , GY JI 

y 

Y ( I ) 

Gyy 

GYY 

a 

ALPHA ( I, J) 

h 

H 

B 

BETA (I , J) 

2hi 

HH  (1,1) 

8x 

BETAX (I , J) 

ki 

PK (1,1) 

By 

BETAY (I, J) 

m 

NPTS 

A9  (cylinder) 

DELTHE 

n 

NSPTS 

(surface) 

DELX 
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Text 

Computer 

Program 

Text 

Computer 

Program 

511 

PHASE1 (1) 

U 

AMU 

612 

PHASE 1 ( 2 ) 

yk  (v) 

AMU  (K) 

*21 

PHASE 2 (1) 

yk(4v) 

AMU 4 (K 

622 

PHASE2 (2) 

V 

ANU 

C 

X(J) 

4v 

ANU4 

n 

Y(J) 

TT 

PI 

IV.  DISCUSSION  AND  RESULTS 


— 

I 

A.  SELECTION  OF  COMPUTER  PROGRAM  PARAMETERS 

The  computer  program  was  developed  on  the  IBM-360 
computer  using  FORTRAN  H.  All  subroutines  were  compiled  on 
DATA  CELL  to  minimize  compile  time  and  core  size.  In  order 
to  maximize  accuracy  while  minimizing  computational  time, 
two  key  computer  parameters  were  evaluated  to  determine 
optimum  values:  cylinder  nodal  points  and  free  surface  model 
points . 

Using  the  first-order  solution  only,  the  number  of  nodal 
points  on  the  cylinder  surface  was  varied  from  12  to  60  and 
compared  with  the  results  determined  by  Garrison  [Ref.  2]. 
The  minimum  number  of  cylinder  surface  nodal  points  that 
provided  accuracy  within  one  percent  of  those  obtained  by 
Garrison  using  60  points,  was  24.  Accordingly,  the  good 
correlation  from  60  points  down  to  24  points  led  to  the 
selection  of  24  nodal  points  on  the  cylinder  surface. 

Since  the  free  surface  integral  was  to  be  carried  out 
from  -oo  to  +°o , the  next  step  was  to  establish  the  finite 
interval  of  convergence  on  the  surface  as  well  as  the  sub- 
division size.  After  the  outer  limits  of  the  free  surface 
integral  were  determined,  the  subdivision  size  was  varied 
from  4 to  126  subdivisions  per  second-order  wave  length 
holding  the  total  interval  constant.  Above  64  subdivisions 
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the  results  did  not  vary  more  than  one  percent,  leading  to 
the  selection  of  64  subdivisions  per  second-order  wave  length 
on  the  free  surface  for  purposes  of  making  calculations. 

Additionally,  the  total  interval,  2a,  was  varied  from 
one  to  eight  wave  lengths  in  both  the  positive  and  negative  x 
directions.  The  second-order  force  coefficients  and  their 
respective  phase  shift  angles  converged  to  a small  value 
varying  periodically.  Convergence  to  this  limit  cycle 
occurred  between  the  first  and  second  wave  length  from 
the  center  of  the  cylinder.  Therefore,  a total  interval  of 
from  -2  to  +2  second-order  wave  lengths  was  chosen  as  adequate 
for  generation  of  the  numerical  results. 

Although  the  values  of  the  second-order  force  coefficients 
and  their  respective  phase  shift  angles  tended  to  converge 
with  increasing  the  limits  of  the  free  surface  integration, 
there  remained,  in  general,  the  small  limit  cycle  as  noted 
above.  The  amplitude  of  the  cycle  was  a function  of  the 
parameter,  a,  and  was  of  significant  magnitude  only  for 
values  of  a less  than  0.25.  The  effects  of  a on  the  limit 
cycle  for  one  of  the  second-order  parameters,  horizontal  force 
coefficient,  is  demonstrated  in  Fig.  (2).  This  represents  a 
plot  of  percent  variation  of  the  force  coefficient  from  the 
mean  versus  the  surface  interval  size,  A/(L/2) , corresponding 
to  values  of  a of  0.25,  0.20  and  0.15.  Figure  (2)  is 
represenLative  of  the  behavior  of  all  second-order  force 
coefficients  ind  phase  shift  angles  and  demonstrates  the 
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SURFACE  INTERVAL  X/(L/2) 

Figure  2.  SECOND- ORDER  HORIZONTAL  WAVE  FORCE 
COEFFICIENT  VERSUS  SURFACE  INTERVAL 


■ 


f 

increase  in  the  amplitude  of  the  variation  with  a decrease 
in  a.  The  amplitude  of  the  periodic  variation  ranges  from 
up  to  fifty  percent  at  a wave  number  of  0.15  and  depth  of 
d = 1.4,  to  less  than  two  percent  for  all  values  of  wave 
number  above  0.25. 

B.  RANGE  OF  APPLICABILITY 

In  a complex  computer  program  of  the  type  developed  to 
carry  out  the  present  numerical  scheme,  certain  limits  with 
respect  to  numerical  stability,  etc.  show  up.  These  computing 
limits  arise  for  various  reasons  such  as  overflows  caused  by 
computing  hyperbolic  functions  of  very  large  numbers, 
numerical  convergence  instabilities,  etc.  While  no  attempt 
was  made  to  investigate  each  of  these  numerical  problems,  a 
list  of  the  limitations  of  the  particular  program  developed 
in  this  work  to  carry  out  the  computations  are  as  follows: 

a — dependent  upon  acceptable  error,  but  0.25 
or  less 

a — deep  water  condition  reached 
h - 3 

h — deep  water  condition  reached  for  a > 0.25, 
h = 20  for  a < 0.25 

d — from  1.4  at  h equal  to  or  less  than  5 down 
to  2 . 5 at  h = 20 

SMIN  — 0.12  when  cylinder  depth  is  other  than 
near  the  free  surface 

SMIN  — 0.30  when  cylinder  depth  is  near  the 
free  surface 


L. 


minimum 

maximum 

minimum 

maximum 

minimum 

minimum 

maximum 
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C.  RESULTS 


A representative  set  of  results  for  the  first-order  and 
second-order  force  coefficients  have  been  generated  for  a 
water  depth  of  h = 5.0,  submergence  depths  of  d = 1.5,  2.0, 

2.5,  and  3.0,  with  a ranging  from  0.15  to  1.2.  All  of  the 
numerical  results  presented  herein  were  based  on  24  sub- 
divisions on  the  circular  cylinder  and  64  subdivisions  per 
second-order  wave  length  on  the  free  surface  integration. 

The  surface  integral  was  carried  out  from  x = -L  to  +L  since 
this  was  found  to  be  adequate  for  convergence.  Since  the 
second-order  wave  length  is  half  the  first-order  wave  length, 

L,  the  total  interval,  is  equal  to  four  second-order  wave 
lengths . 

The  first-order  horizontal  and  vertical  force  coefficients 
are  presented  in  Figs.  (3)  and  (4),  respectively.  It  may  be 
noted  that,  in  general,  the  forces  decrease  with  depth  of 
submergence  according  to  expectation  since  the  wave  action 
dies  out  with  depth. 

The  second-order  horizontal  and  vertical  force  coefficients 
are  presented  in  Figs.  (5)  and  (6),  respectively.  Generally, 
the  results  show  the  second-order  effect  to  be  relatively  more 
important  as  the  cylinder  approaches  the  free  surface.  It 
might  be  suspected  from  this  that  the  second-order  contribution 
would  be  much  more  significant  in  the  case  of  surface  piercing 
or  floating  bodies. 
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The  horizontal  and  vertical  steady-state  force  coefficients 
are  shown  in  Figs.  (7)  and  (8),  respectively.  These  results 
show  that  the  horizontal  steady-state  force  coefficient  is  very 
small  in  general.  The  horizontal  force  can  be  shown  to  be 
proportional  to  the  momentum  flux  of  the  reflected  wave  and 
the  reflected  wave  is  small.  In  fact,  in  the  case  of  infinite 
water  depth.  Dean  [Ref.  1]  showed  that  the  reflection  was 
exactly  zero,  and,  accordingly,  the  steady-state  horizontal  force 
is  likewise  zero. 

It  may  be  noted  that  the  steady-state  vertical  force  is 
positive.  This  results  from  the  fact  that  the  velocities  are 
largest  on  the  top  of  the  cylinder  and,  accordingly,  the  average 
pressure  is  reduced  on  the  top  of  the  cylinder. 

The  phase  shift  angles  of  the  first-order  and  second-order 
forces  are  shown  in  Figs.  (9)  - (12)  . 

To  demonstrate  the  effect  that  the  second-order  terms  have 
on  the  forces  on  the  cylinder  and  their  respective  phase  shift 
angles,  a comparison  was  made  between  the  first-  and  second- 
order  results.  Figures  (13)  - (18)  are  plots  of  the  horizontal 
and  vertical  dimensionless  forces  on  the  cylinder  versus  time 
over  one  complete  wave  cycle  for  both  the  first-order  solution 
alone  and  the  combined  first-order  and  second-order  effects. 
Additionally,  a plot  of  the  incident  wave  is  included  to 
demonstrate  the  phase  shift  involved  (the  amplitude  of  the 
incident  wave  has  no  significance).  A mean  water  depth,  h,  of 
5.0,  a cylinder  depth,  d,  of  1.5  and  a wave  height,  H,  of 
0.5  were  used,  with  the  wave  number,  a,  assigned  three  values, 
0.25,  0.5  and  1.0. 
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WAVE  NUMBER,  a = 2 7ra/L 

Figure  6.  SECOND-ORDER  VERTICAL  WAVE  FORCE  COEFFICIENT ; 
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WAVE  NUMBER,  a=  2ttci/L 

Figure  9.  FIRST-ORDER  HORIZONTAL  PHASE  SHIFT  ANGLE 
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First-Order 
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TIME,  t (rad) 

Figure  13.  HORIZONTAL  WAVE  FORCE:  a = 0.25,  h = 5.0,  d = 1.5,  H = 0.5 


Second-Order 


igure  14.  VERTICAL  WAVE  FORCE = a =0.25,  h = 5.0,  d = 1.5,  H = 0.5 
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0 1.0  2.0  3.0 

TIME 

Figure  17.  HORIZONTAL  WAVE  FORCE : 
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V.  CONCLUSIONS 


A computer  program  has  been  developed  to  carry  out  the 
numerical  solution  of  the  second-order  wave  — cylinder 
interaction  problem.  The  first-order,  second-order,  and 
steady-state  force  coefficients  were  determined  for  the 
submerged  horizontal  cylinder. 

An  approximate  method  for  dealing  with  the  second-order, 
nonhomogeneous  free  surface  boundary  condition  was  developed 
which  appears  to  converge  except  at  small  values  of  a,  i.e., 
at  very  large  wave  lengths. 
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APPENDIX  A 

COMPUTER  PROGRAM.  LISTING 

THIS  PROGRAM  CALCULATES  WAVE  FORCES  AND  THEIR  PHASE 
SHIFT  ANGLES  FOR  WAVE  INTERACTION  WITH  A SUdMERGEO 
HORIZONTAL  CYLINDER 

A = 2*P  I*CYL INOER  RADIUS/WAVE  LENGTH  = SIGMA**2* 
CYLINDER  RADIUS/ACCELERATIuN  OF  GRAVITY 
H = MEAN  WATER  DEPTH/CYLINDER  RADIUS 
D = CYLINDER  DEPTH/CYLINDER  RADIUS 

NPTS  = NUMBER  OF  CYLINDER  SURFACE  ELEMENTS  AND  NODAL 
POINTS  FOR  NUMERICAL  EVALUATION 

NSPTS  = NUMBER  OF  FREE  SURFACE  ELEMENTS  AND  NODAL 
POINTS  FOR  NUMERICAL  EVALUATION 


COMPLEX  GI  J,GNI  J,GNJ I ,GXI  J,  GXJ I , GY  I J , GY J I , G I J EX T , GYY , 
1DET,SUM1,SUM2,S'JM3,SUM4,U1I  S,U1ISX,U1ISY,U11SYY, 
2U1SS»U1SSX,U1SSY  tUlSSYY 
COMPLEX  ALPHA(24,24),B£TA(24,24),BETAX(24,24), 

1BETAY  ( 24,  24)  ,Hr)(  24,1)  ,F  < 24,  1)  , PK  ( 24,  1)  , Fl(  24, 1)  , 
2FSI500)  {Ul(24),UIX(2<t),UIY(24)  ,U2(24)  , J<:SC  1 ( 24)  , 
3U2SCN1 ( 24 J ,U2SCO<  24) 

DIMENSION  XI 24) , Y ( 24) ,ANX( 24) , ANY! 24) ,CHY( 25) ,CHY2<  25) 
1 ,SHY(  25) , SHY 2(25), COE FG1200), COE FG2( 200 ) , 

2COSAMU  ( 200,2  5)  , COS AM2(  200 ,25) .SINAMUI  200  ,25), 

3 S I NAM  2 ( 20  0,25) , AMU ( 20  0) , AMU41200) ,SH2Y( 24) , CH2Y( 24) 
DIMENSION  Cl (2) ,C2(2) ,C3(2) ,PHASE1(2) ,PHASE2( 2) 

COMMON/ CPX/HH,?K 
CCMMON/VAR/X, Y, A NX, ANY 

COMMCN/GSHY/CHY , CHY2 ,SHY, SHY2,SH2Y,CH2Y 
CCMMON/GMU/COSAMU ,C0SAM2 ,SI NAMU, SINAM2 , AMU , AMU4 , COEFG, 
1C0EFG2 

COMMON/ CON  ST /H,D, DELTHE ,SMI N, NPTS, NSPTS 
COM MON/CP  I / P I 
EQUIVALENCE  ( Hh ( 1 ) , F ( 1 ) ) 

EQUIVALENCE  ( PK ( 1 ) , F L ( L ) ) 

PI=3. 141592 


READ  INPUT  DATA  AND  CALCULATE  REQUIRED  REPEATING 
DATA  AND  ARRAYS 


CALL  GEODAT  ( A,  A2 , ANU,ANU4»  SH2AH, SH2AH2 , SHAH, SHAH2, 
1CHAH,CHAH2 , AO, A A , 3B , CC , DD , A02 , A A2 , BB2 , C C 2 , C D2 ) 

DO  100  1 = 1, NPTS 
DO  100  J=1,I 
XV=X(  I > 

Y V=Y ( [ ) 

XC  = X ( J) 

Y C = Y ( J ) 

ANXI  = ANX(  I ) 

ANY  I = AN Y ( I ) 

ANXJ  = ANX ( J ) 

ANY J = ANY( J ) 

IF(ABS(X(  D-Xl  J)  ) .LT.SMIN)  GO  TO  50 
SHYI  = SHY ( I ) 

C HY I =CH  Y ( I ) 

SHY J=SHY( J ) 

CHY J = CH Y ( J ) 

CALL  GREENS  I A , ANU , S H2 AH, SHAH, CHAH, CCS AMU , S INAMU , AMU , 
1CQEFG,SHYI,CHYI,SHYJ,CHYJ,1,J,XV,YV,XC,YC,ANXI,ANYI, 
2ANXJ,ANYJ,GIJ,GNIJ,GNJI,GXIJ,GXJ I,Gy: J,GYJl , GYY ) 

GC  TO  90 
50  CONTINUE 

CALL  GREEN  ( A , ANU , SH2 AH , SHA H , CH AH , AO , A A , B8 , CC , DD , I , J , 
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r 


1XV,YV,XC, YC, ANXI , ANY I,GIJ,GNIJ,GXIJ,3YIJ,GYY) 

IF(I.EQ.J)  GN J I =GN  I J 
IF  (I.EQ.J)  GX J I =GX I J 
IF  (I.EQ.J)  GY  J I =GY I J 
I F ( I.EQ.J)  GO  TO  90 

CALL  GREEN  ( A , ANU , SH2 AH , SHA H.CHAH , AO , A A , BB , CC , DD , J , 1 , 
' 1JXC,  YC  ,XV,  YV  , ANX  J , ANYJ.GI  JEXT,  GNJ  I , GX  J I , G Y J I , GY  Y ) 

90  CONTINUE 


CALCULATE  THE  >1  RST-ORDER  ALPHA  AND  BETA  MATRICES 


ALPHA { I *J)=( l./PI ) *GNI J*DELTHE 
ALPHA( J f I ) = (!./ PI ) *GN  J I *DEL  THE 
8ETA( I , J)  = ( I./ { 2 .*PI ) )*GI J* DEL  THE 
B ETA ( J , I) -BETA ( I , J) 

BETAXd  , J)  = ( l./(2.*PI  ) )*GXI  J*DELTHE 
BET AX (J  , I )=(  l ./(  2 . *P  I ) 5*3XJI*OELTHE 
B ETAY ( I,J)=(l./(2.*PI)  )*GYI J*DEITHE 
B E TAY ( J » I ) = ( 1./ ( 2 . *PI ) ) *GY  J I^OELTHE 
100  CONTINUE 

00  120  1 = 1 j N PTS 

ALPHA ( I , I) =ALPHA( I , I ) +CMPLX ( 1.0,0 .0) 

BETAX(  I , I)=BETAX ( I , I) +ANX ( I )*CMPLX (0.5, 0 .0) 
BET  AY  (I  • IJ=BETAYC  I • I )+■  ANY  ( I )#CMPLX(0.5,  0.0) 
120  CONTINUE 


GENERATION  OF  THE  FIRST-ORDER  DISTRIBUTION  FUNCTION, 
F ( I , 1 ) , BY  INVERSION  OF  ALPHA!  I » J ) *F( i , 1 } = HH (1*1) 


CALL  COMAT  l 2A- , 1 , ALPHA, HH,DET,  INDICA) 


THE  HH (1,1)  MATRIX  IS  REPLACED  BY  THE  DISTRIBUTION 
FUNCTION, F< I ,1 ) 

COMPUTATION  OF  THE  FIRST-ORDER  SCATTERING  PuTENT I AL 
FUNCTION,  U1SC(I),  AND  ITS  X AND  Y PARTIAL  DERIVATIVES 
U 1 SCX (II  AND  JiSCYlI),  AND  COMBINATION  WITH  THE 
INCIDENT  POTENTIAL  COUNTERPARTS  TO  COMPUTE  THE  TOTAL 
POTENTIAL  FUNCTION  AND  ITS  X AND  Y DERIVATIVES 


DO  155  1 =1 , NPT  S 
SUM 1=( 0.0, 0.0) 

SUM2=<  0 .0,0.0) 

SUM3=(0. 0,0.0) 

DO  150  J= 1 » N PT S 
SUM1=SUMUF(  J,  1)  *BETA<  I , J) 

SUM2=  SUM2+F ( J,1)*BETAX(I,J) 

SUM3=SUM3+F( J, 1 ) *B  ET  AY ( I , J ) 

150  CONTINJE 

Ul( I )=SUM1-(CHY ( I ) / ( A *CH AH ) )*CEXP(CMPLX(0.0,A*X(I))) 
U1X(I)=SU<M2-CMPLX(  0.0,1.0)*CHY  ( I ) *CEXP(  CMPLX(  0.0, 

1 A*X ( I ) ) J/CHAH 

U 1 Y ( I ) = SUM3-SHY ( I ) *CEXP(CMPLX( 0.0  ,A*X( I ) ) ) /CHAH 
155  CONTINUE 


EVALUATION  OF  THE  FREE  SURFACE  PRESSURE  DISTRIBUTION 
SOURCE  STRENGTH  FUNCTION,  FS(L) 


DELX=0.015625*PI/A 
kvRITE  ( 6 , 1 BO  ) DE  LX 
130  FORMAT  ( 5X/ / 5X , • D E LX  = • , F 1 0 . 5 / / / ) 
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SP=0.0 

CST=.5*DELTHE/PI 
I COUNT  = 1 

DU  250  L=1 ,NSPTS 
SUM  1= ( 0 .0,0.0) 

SUM2=( 0 .0,0.0) 

SUM3=(0. 0,0.0) 

SUM<t=(0  .0,0.0) 

DO  245  1=1 , N PTS 
XV=SP 
V V=Q . 0 
XC=X( I ) 

YC=Y( I ) 

A NX  I =0 . 0 
ANY  I = 1 . 0 
ANXJ=ANX<  I ) 

ANY  J=AN Y ( I ) 

SHY  1=  SHY ( 2 5 ) 

CHYI=CHY( 25) 

SHY J=SHY( I ) 

C HY  J = CH Y ( I ) 

IF  IABSIXV-XC) .LE.SMIN)  GO  TO  240 

CALL  GREENS  ( A , ANU , SH2AH, SHAH, CHAH.COSAMU, S INAMU , AMU, 
1C0EFG,SHYI , CHY i , SHY  J , CHY  J , 2 5 , I , X V , YV , XC , YC , ANX I , A NY  I , 
2ANXJ, ANYJ,GI J,GNI J,GNJI , GX I J , GX J I , GY  I J , G Y J I , GYY ) 

GO  TO  241 

240  CONTINUE 

CALL  GREEN  ( A , ANU  , SH2 AH , SHA H,CH AH , AO , AA , BB , CC , 00 , 2 5 , I , 
1XV,YV,XC,YC, ANX I , ANY  I , G I J , GN 1 J , GX I J , GY  I J , GY Y ) 

241  CONTINUE 

SUM 1 = SUM 1 +CST#GIJ*F(  1,1) 

SUM2=SUM2*CST*GX I J*F  ( 1,1) 

SUM3=SUM3«-CST*GY  I J*=F  I 1,1) 

SUM4=SUM4*CST*GYY*Fi  1,1) 

245  CONTINUE 

U1IS=(-1./A)*CEXP(CMPLXIO.O,A*SP)  ) 

U1  ISX=CMPLXt 0. 0, -1 .0) *CEXP( CMPLX10.0 , A*SP)  ) 
U1ISY=-TANH< A*H) *CEXP(CMPLX (0.0 ,A*SP) ) 

U II  SYY  =-A*CEXP(CMPLX(0.0,A*SP)  ) 

U1SS=SUM1 
U 1 SSX=  SUM2 
U 1 SSY=  SUM3 
U1SSYY=SUM4 

FS(L)  = ( 2. * A/ (3.* ANU ) )*<  J1 I S*U1 SSYY+U1 ISYY^UISS+UISS* 
1U1SSYY-6.*U1  ISY*UISSY-3.*U1  SSY*UlSSY-4.  *Ul  I SX-U1  SSX 
2-2.*UlSSX*UlSSX) 

IF  ( ICOUNT .EQ.O ) GO  TO  248 
SP=FL0AT(L+l)*DELX/2. 

I CCUNT=0 
GO  TO  250 
248  SP=-SP 

I COUNT  = 1 
250  CONTINUE 


CALCULATION  OF  THE  PARTICULAR  SOLUTION  PORTION  OF  THE 
SECOND-ORDER  SCATTERING  POTENTIAL  AND  ITS  NORMAL 
DERIVATIVE,  U2SCKI)  AND  U2SCNKI) 


CST=.5*DELX/PI 
DO  350  1=1 ,NPTS 
SUM  1= ( 0.0,0.01 
SUM2=(0. 0,0.0) 
SP=0.0 
I COUNT  = 1 

DO  345  L=1 »NSPTS 
XV  = X(  I ) 
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1 


YV  = Y( I ) 

XS  = SP 
Y S = 0. 0 
ANX 1 = ANX ( I ) 

ANY  I = AN Y ( I ) 

ANX J=0 . 0 
ANY  J= 1 . 0 
SHY  1 = SHY ( I ) 

CHY I = CH Y ( I ) 

SHYJ=SHY<  25) 

CHYJ=CHY{ 25) 

IF  (ABS(XV-XS) .LE.SMIN)  GO  TO  340 

CALL  GREENS  ( A2 , ANU4 , SH2AH2 , SHAH2 ,CHAH2 , C0SAM2 , S I NAM2, 
1AMU4,C3EFG2*  SHY  I , CHY I ,SHY J,CHY J , I , 25 , XV , YV , XS , YS ♦ 
2ANXI,ANYI, ANXJ , ANY J,GI J ,GNI J,GNJI,GXI J,GXJI ,GYIJ,GYJI, 
3GYY) 

GO  TO  341 

340  CONTINUE 

CALL  GREEN  ( A2 , ANU4 , SH2 AH2 , SHAH2 , CHAH2 , A 02 , AA2 , BB2 , CC2 
1 , DO 2, I , 2 5, XV, YV, XS,YS,ANXI , ANY  I , G I J , SNI J , GX  I J , GY  I J , 
2GYY ) 

341  CONTINUE 

SUM1  = SUM1*CST*GI J*FS(  L) 

SUH2«SUH2*CST*GNI J*FS(  L) 

IF  ( ICOUNT.EQ.O)  GO  TO  290 
SP=FLOAT( L+l )*DELX/2 . 

I COUNT=  0 
GO  TO  345 
290  SP=-SP 

I COUNT=  1 
345  CONTINUE 

U2SC1( I ) =SUM  I 
U2SCN1 ( I ) = SUM2 
350  CONTINUE 


CALCULATION  OF  THE  HOMOGENEOUS  SOLUTION  PORTION 
OF  THE  SECOND-ORDER  SCATTERING  POTENTIAL,  U2SC0(I) 


DO  500  1 = 1 , NPT  S 
DG  500  J=1,I 
XV=X( I ) 

Y V = Y l I ) 

XC=X( J ) 

YC=Y(J ) 

ANX I = ANX ( I ) 

ANYI  = AN Y ( I ) 

ANXJ=ANX( J ) 

ANY J = ANY ( J ) 

IFIABSIXI I )-X( J) ) .LT.SMIN)  GO  TO  450 
SHY  I = SH  Y ( I ) 

CHY I = CH Y ( I ) 

SHY  J=  SH  Y ( J ) 

CHY  J=CH Y ( J ) 

CALL  GREENS  ( A2 , ANU4, SH2AH2 , SHAH2 ,CHAH2 , C0SAM2, S INAM2, 
1AMU4,C0EFG2, SHY I,CHYI,SHYJ,CHYJ,I,J,XV,YV,XC,YC,ANXI, 

2 ANY  I , ANXJ, ANY J, GI J ,GNI J.GNJ I , GX I J , GX J I , G Y I J , GY J [ , G YY  ) 

GC  TO  490 
450  CONTINUE 

CALL  GREEN  < A2  , ANU4 , S H2 AH2 , SHA H2 , CHAH2 , A 02 , AA 2 , 3 B 2 , CC2 
1 »0D2, I,J,XV,YV,XC,YC, ANX I , ANY  I , G I J , GN I J , GX I J, GY  I J , 

2GYY ) 

IF(I.EQ.J)  GN J I =GN  I J 
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I F ( I.EQ.J)  GO  TO  4-90 


CALL  GREEN  l A2 , ANU4 , S HZ AH2 , SHAH 2 , CHAH2 , A02  , AA 2 , BB 2 , CC2 
1 ,CD2, J, I ,XC,YC ,XV, YV ,ANXJ, ANY J,GI JEXT,GNJ I ,GX JI ,GYJ I , 
2GYY) 

490  CONTINUE 

ALPHA! I , J) =( l./PI ) *GY I J *D  EL  THE 
ALPHA!  J , I } = ( I ,/P  I )*GYJI*DELTHE 
BETA(I,J)=!1./(2.*PI))*GIJ*DEITHE 
BETA! J,  I ) = BE  T A ( 1 ,J) 

500  CONTINUE 

OG  510  1=1 ,NPTS 

PK( 1 , 1) = 2*( PK( 1 , 1J-U2SCN1 ( I ) ) 

510  CONTINUE 

DO  520  1=1 ,NPTS 

ALPHA! I , I) = ALPHA! 1,1) +CMPLX ( 1 .0, 0.0) 

520  CONTINUE 

CALL  COMAT  ( 24 , 1 , AL PH A , PK , D E T , I ND I C A ) 

DO  540  1=1 , NPTS 
S UM= ( 0.0, 0.0) 

DO  530  J = 1 ,NPTS 
SUM=SUM+F 1!J,1)*BETA(I,J) 

530  CONTINUE 

U2SC0 ( I ) =SUM 
540  CONTINUE 


EVALUATION  OF  THE  TOTAL  SECOND-ORDER  POTENTIAL,  U 2 ( I ) 


DO  550  1=1 , NPTS 

U2( I )=U2$C1( I ) +U2SC0!  I )-(CHY(  I ) / (2.* A* ( SHAH**4) ) ) * 
1CEXP(CMPLX(C.0,2.*A*X(  I))) 

550  CONTINUE 

WRITE  ( 6,560) 

560  FORMAT  ( 5 X/ / 9X , 1 I ' , 1 5 X , 1 U 1 ( I ) ' , 24X  , • U 1 X ! I ) ' , 2 5X  , 
1'U1Y(I  ) ' , 24X , ' U 2 ( I )'//) 

00  530  1=1, NPTS 

WRITE  1 6,570)  I , U 1 l I ) , U 1 X ( I ) , J 1 Y l I ) , J 2 ( I ) 

570  FORMAT  ( 5X , I 5 , 4 ! F 1 6 . 6 , F 14 . 6 ) ) 

580  CONTINUE 


EVALUATION  OF  THE  FIRST-ORDER,  SECOND-ORDER  PERIODIC, 
AND  SECOND-ORDER  STEADY  STATE  FORCE  COEFFICIENTS  AND 
THE  PERIODIC  FORCE  PHASE  SHIFT  ANGLES 


SUM1=(0. 0,0.0) 

SUM2= ( 0 .0,0.0) 

SUM3=(0. 0,0.0) 

SUM4= ( 0 .0,0.0) 

SUM5=0 . 0 

SUM6=0. 0 

DO  720  1 = 1 , NPTS 

SUM1  = SUM1 +U1  ( I ) * A NX!  I ) vQELTHE 

SUM2=SUM2*-Ul  1 I) *ANY( I > *D£L  THE 

S JM3=SUM3+ l ( 6.  * ANU*ANU*U2  !I)/A)-U1X(I)*U1XI  I1-U1Y1I) 
1 *U1 Y ( I ) )*ANX< I ) *DELTHE 

SUM4=SUM4+< ( 6.  *ANU*ANU*U2 ( I)/A)-ULX(I)*U1X<  I ) — U l Y ( I ) 
1 *U 1 Y < I ) ) *ANY ( I ) *Q  EL  T HE 

SUM5  = SUM5M CAES ! U 1 X ( I)*U1X!I))*CA8S!U1YI I )*U1Y( I ) ) 
H-ANU*ANU/A**2-1 . 0 ) *ANX ( I ) *D EL  THE 
SUM6=SUM6* l CABS ( U 1 X ( I ) *U1X ( I))+CA6$(U1Y(  I ) *U1Y{ I ) ) 

1 «-  ANU*ANU/A**2~1  . 0 ) * AN  Y ( I ) *D  EL  T HE 
720  CONTINUE 

Ci (1 )=CABS( SUM1 * A ) 

C 1 < 2) = CABS ( SUM2  * A ) 

C2ll)=CABS(SUM3*A*A/<4.*ANU)  ) 
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C2(2)=CABS(SUM4*A*A/(4.*ANU) ) 

C 3 ( 1) =SUM5*A*A/( 4. *ANU) 

C3(2)=SUM6*A*A/(4.*ANU) 

PHASE  1 ( 1 ) = ATAN2  ( A I MAG  ( SUM1*A)  ,REAL(SUM1*A)) 

PhASEl ( 2)=ATAN2l  A I MAG ( SUM2*A) ,REAL(SUM2*A)) 

PHASE  2 1 1 ) = AT  AN2 ( A I MAG ( SUM3* A*A/ ( 4 . *ANU ) ) , RE AL ( SUM  3* A* A 
1/ ( 4.*ANU) ) ) 

PHASE2 I 2)  = AT AN2 ( A I M AG ( SUM4* A *A/ ( 4.*ANU)  ) , RE  AL ( SUM4*A*A 
1/ (4.*ANU) ) ) 

WRITE  (6,730)  C 1 ( 1) ,C1(2)  ,C2(  1)  »C2(2) .03(1)  ,C3(2), 

1 PHASE  1 ( 1 ) , PHASE  1 (2>,PHASE2(1),PHASE2(2) 

730  FORMAT  l 5X / / / 5 X , • C 1 ( 1 ) = ' , F8 . 5 , 5X  , • Cl  ( 2 ) = ' , F 8 . 5/ / 5X  , 

1 • C2( 1) =• , F8. 5, 5X, • C2( 2)=’ , F8.5//5X,  *C3(  1 ) = ' ,F8.5, 5X, 

2'  C3(2)  ='  .F8.5//5X,  • PHASE  1 ( 1 ) = ' , F8.5.5X,  ' PHASEK2)  =•  , 
3F8.5//5X, 'PHASE 2(1)='  ,F8. 5, 5X,  ' PHASE2<2)=' ,F8.5//) 

STOP 

END 


THIS  SUBROUTINE  READS  THE  INPUT  GEOMETRICAL  DATA  AND 
CALCULATES  THE  FIRST  200  RCCTS  OF  AMU*T AN( AMU*H)  + 

ANU  = 0 AND  AMJ4*TAN( AMU4*H  ) ANU4  =0.  IT  ALSO 
GENERATES  ARRAYS  OF  CERTAIN  FUNCTIONS  AND  COEFFICIENTS 
USED  IN  GREEN  AND  GREENS  AS  WELL  AS  THE  Hh  AND  PK 
MATRICES 


SUBROUTINE  GEODAT  ( A , A2 , ANU , ANU4, SH2AH , SH2AH2 , SHAH , 
lSHAH2,CHAh,CHAH2 , A0,AA,B6,CC,DD, A02,Aa2 , BB2 ,CC2, DD2) 

COMPLEX  HH ( 24, 1 ) , PK ( 24, 1 ) 

DIMENSION  X( 24) , Y( 24) , ANX( 24) , ANY ( 24) ,CHY(25),CHY2(25) 
1 , SHY (25) , SHY2( 25) ,CCEFG( 200) ,COEFG2( 200 ) , AMU( 200) , 

2 A MU4( 20  0) .COSAMJ ( 20C ,25) ,COSAM2 (200,25), SINAMU(  200,25) 
3, SIN A M2 (200, 25 ) »SH2Y{ 24) ,CH2Y( 24) ,XX(  24) 

CCMMQN/CPX/hh, p< 

COMMON /VAR /X, Y, ANX.ANY 

COMMON/ GSHY/C HY ,CHY 2, SHY , SHY  2 , SH2Y , CH2Y 
C0MM0N/C-MU/C0SAMU,CUSAM2  , SI  NAMU,  SINAM2  , AMU  , AMU4 , CCE  FG, 
1C0EFG2 

COMMON/C 0NST/H,D,D6LTHE,SM I N.NPTS, NS PTS 
COMMON /CP  I /P  I 

RE AD (5, 5)  A,H,D, SM  IN, NPTS , NSPT S 
5 FORMAT  ( 4F 10.5,214) 

ANU=A*TANH( A*H) 

SH2AH=S INH( 2 . * A*  H ) 

CHAH=COSH( A*H ) 

SHAH=  S I N H(  A*H ) 

AO  = T AN  H( A*h) 

AA=  1 . /CHAH**2 
E6=-H*SHAH/(CHAH**2) 

CC =-H*H*<  1 .-SHAH* *2)/  ( 3.*CHAH**3) 

DD  = H**3*SHAH*(  2.*CHAH**2*3.*(  1 . - SHAH**  2 ) ) / ( CHAH**4* 12) 

0ELTHE=2.*PI/NPTS 

THE  TA  = DEL  THE/2 . 

DO  15  1=1, NPTS 
X (I)=COS(THETA) 

Y ( I ) = Si THET  A ) -D 
A NX  ( I ) = X ( I ) 

ANY ( I ) = Y ( I )*D 

SHY ( I ) = S INH( A* ( H*Y ( I ) ) ) 

C HY  < I )=COSH( A* ( H + Y ( I ) ) ) 

SH2Y ( I)=SINH(2.*A*(H+Y(  I ) ) ) 

C H2Y ( I ) = COSH<  2.* A* (h*Y( I ) ) ) 

HH(  I , 1 ) = ( 2 ,/CHAH)  *CMPLX(  ANY  (I)*SHY(I),ANX(I  )*CHY(  I ) )* 
1CEXP( CMPLX ( 0.0 , A*X l I ) ) ) 

PK(  1, 1 ) =(  ANY(  I ) *SH2Y  ( I ) *CMPLX  ( 0 . 0 , l . 0)  * A NX  ( I)  *Cri2Y(  I ) ) 
1*CEXP(CMPLX(0»C,2 ,0*A*X( I ) ) ) /SHAH**4 
THETA=THETAfDELTHE 
15  CONTINUE 

SHY ( 2 5 ) = SH AH 
CHY<  25) =CHAH 
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6=ANU*H 
DO  50  K= 1 » 200 
XX( 1 >=P1*K 
DO  25  1=1,20 
I 1 = 1 

Y Y=XX ( I ) 

XX ( 1*1) =XX ( 1 ) - A T AN2 ( 3 , YY ) 

1 F ( A3S( ( XX( I +1 ) -XXI I) ) /XX( I +1 ) ) .LT. .0001 ) GO  TO  30 
25  CONTINUE 

IF  ( II.GE.20)  GO  TO  27 
GO  TO  30 

27  WRITE  (6,28) 

28  FORMAT  (5X,'AMU  ROOT  DOES  NOT  CONVERGE'//) 

30  CONTINUE 

AMU(K) =XX  (II  ) / H 

COEFG (K)  = 2.*PI*(  AMUU)  *AMU<  K) *ANU*ANU) / ( ANU*AMU(  K) -H* 
1 AMU(K) *( AMU(K) *A,MU ( K) +ANU-ANU ) ) 

NPTS1=NPTS+1 
DO  AO  I = 1 , NP  T S 1 
IF  ( I .EC.N°T51 ) GO  TO  35 
COSAMUl K,I )=COS(AMU(K)*(H*Y(I) ) ) 

S INAMU(  K,  I)=SIN(AMU(X)*t(H  + Y(  I ) ) ) 

GO  TO  AC 

35  COSAMU(K,I )=COS(AMU(K)*H) 

S INAMU( K,l)=SiN(AMU(K)*H) 

AO  CONTINUE 
50  CONTINUE 

ANUA=A. *ANU 
B2=ANUA*H 
DO  70  K= 1 , 200 
XX( 1 ) = P I *K 
DO  55  1=1, 20 
I 1 = I 

Y Y = XX ( I ) 

XX ( 1 + 1 ) =XX( 1 )-ATAN2( 82  ,YY  ) 

IF( A8S(  (XX( 1+1 ) -XX(  I ) ) /XX( 1*1) ) .LT. .0001  ) GO  TO  60 
55  CONTINUE 

IF  ( II.GE.20)  GO  TO  57 
GO  TO  6C 

57  WRITE  (6,58) 

58  FORMAT  (5X,'AMUA  ROOT  DOES  NOT  CONVERGE'//) 

60  CONTINUE 

AMUA(K) =XX( I I ) / H 

CCEFG2 ( K ) = 2.  *P I * ( AMUA ( K ) *AMUA( K) +ANUA*ANUA)  /( ANUA* 

1 AMUA ( K ) -H*AMUA( K ) * ( AMUA ( K ) * AMUA ( K ) +ANUA* ANU A ) ) 
NPTS1=NFTS+1 
DC  65  I =1 , NPTS 1 
IF  (I .EQ.NPTF1 ) GO  TO  68 
C0SAM21 K, I ) = COS ( AMUA ( K ) * ( H+Y < I ) ) ) 

SINAM2 ( K, I )=SI  Nl AMUAI K )* ( H+Yl I ) ) ) 

GO  TO  65 

68  CQSAM2(K,I )=COS(AMUA(K)*H) 

S I NAM2 ( K , I ) = S I N ( AMUA ( K) *H) 

65  CONTINUE 
70  CONTINUE 
XX< 1 )=ANUA 
DO  80  1=1,20 
I I = I 

Y 2 = XX ( I ) 

XX ( I*  l ) =A. *ANU/ T ANH l Y2*H ) 

IF  ( ASS  ( >.  XX(  I +1 ) -XX(  I ) i /XX(  I +1 ) ) .LT. 0.0001 ) GO  TO  85 
30  CONTINUE 

WRITE  (6,3  2) 

82  FORMAT  (5X,'A2  ROOT  DOES  NOT  CONVERGE'//) 

85  CONTINUE 
A 2=  XX ( I I ) 

SH2AH2=SINH( 2.*A2*H) 

CHAH2-C0SH< A2*H) 

SHAH2=  S INH ( A2*H) 

A02  = T ANH  ( A 2 * H ) 

AA2=1./CHAH2**2 
B32=-H*SHAH2/( CHAH2**2) 
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CC2=-H*H*( l.-SHAH2**2) / (3.*CHAH2**3 ) 
DD2=(H**3)*SHAH2*(  2 . *CHAH  2*  * 2 + 3 . * < l.-SHAH2**2))/ 

1 ( ( CHAH2**A) * 12  . ) 

DO  90  I=1,NPTS 

SHY  2(  I ) =S IMH( A2* (H  + Y(  I J ) ) 

C HY 2 ( I ) =COSH( A2* ( H+Y(  I ) ) ) 

90  CONTINUE 

WRITE  (6,95)  A,A2,D,H,ANU,ANU<+,NPTS,SMIN,NSPTS 
95  FORMAT  ( 5X  / / / 3 X , ' A = • , F 1 0 . 5 , , ' A 2 = • , F 1 0 . 5 , A X , • D=  ' , 

1 F 10 . 5 , A X , 1 H= ' , F 1 0 . 5 » 4X , ' ANU=' ,F  10.5,^X , • ANU4=* ,F 10. 5/ 
24X,'NPTS=,,I3,‘+X,'SMIN=',F10.5,4X,'N$PTS  = ,,I3///) 
RETURN 
END 


THIS  SUBROUTINE  CALCULATES  G AND  ITS  DERIVATIVES, 

GX , GY , GN , AND  3YY , 3Y  USE  CF  THE  INTEGRAL  FORM  FOR  THE 
CASE  (XII)  - X(J))  LESS  THAN  SMIN 


SUBROUTINE  GREEN  ( A , A NU , S H2 AH, SH AH, CHAH , AO , AA , SB , CC , DO 
1,I,J,X,Y,XI,ETA, ANX , ANY , G, GN, GX , GY , GYY ) 

COMPLEX  G,  GN 
COMPLEX  GX , GY , G Y Y 

DIMENSION  TEST(  200 ) , TESTT(  100)  , TESTTT(  I 00)  ,SUMOX(  15)  , 

1 SUMOY ( 15) ,NNN( 15),TESTYY(200) 

CCMMON/CONST/H, D, DEL THE, SMI N,NPTS ,NiPTS 
COMMON/CP  I/P  I 

PI (Y, ETA, XX, AMU) =COSH ( AMU* ( Y+H ) )*COSH(AMU*(  ETA+H ) ) * 
1C3S( AMU*XX ) / l COSH ( AMU*H)**2 ) 
P2X(Y,ETA,X,XI,AMU)=-AMU*C0Sri(AMU*(ETA+H  ) ) *SIN(  AMU* 

1 (X-X I ) ) *COSH( AMU* ( Y +H  ) ) / ( COSH ( AMU*H)**2) 
P2Y(Y,ETA,X,XI,AMU  ) =AMU*COSH(  AMU*  (ETA+H ) ) *S  INH(  AMU* 

1 ( Y + H) ) *COS (AMU* (X-X I ) ) / ( C OS H ( A MU*H ) **2 ) 

P2YY ( Y , ETA ,X ,X I , A MU ) = AMU* AMU*COS H ( AM J * ( Y + H ) ) *COSH( AMU* 
1 ( ET  A + H I ) *COS(  AMU*( X-X I ) ) / ( C CSH  ( AMU*H ) ** 2 ) 
w 1 (Y ,ETA,XX, AMJ) = - P 1 ( Y , E T A , X X , AMU ) * ( AMU- A)/ (AMU* 

1 T ANH ( AMU*H)- AN J ) 

Q2X(Y,ETA,X,XI,AMU)=-P2X(Y,ETA,X , XI, AMU ) *(  A MU- A)  / (AMU* 
1 T ANH( AMU*H)-ANJ) 

Q2Y(  Y , ETA  , X,  XI  , AMU  ) =-P2Y(  Y , ETA  , X , XI  , AM'J  ) *(  A MU- A)  / ( AMU* 
1 T ANH ( AMU*H)-ANU ) 

Q2YY(  Y, ETA,X,XI , AMU ) =- P2 Y Y ( Y , £ T A , X , X I , AMU ) * ( AMU- A ) / 

1 ( AMU*TANH( AMU*H ) -ANUJ 

Q1C(Y, ETA, XX )=-COSH( A* (Y+H)  ) *COSH( A* (ET A + H)  )*COS  < A *XX) 
1* A/(COSH( A*H)**2*( ( A*A-ANU*ANU)*H+ANU) ) 

Q2CX l Y , ETA ,X ,XI ) =A*COSH( A*( ETA  + H ) ) *S IN( A* ( X -X  I ) ) *COSH 
1( A*(Y+H) )*A/(C0SH(A*H)**2*(  ( A* A-ANU*ANJ ) *H+ ANU) ) 

Q 20Y( Y , E T A , X , X I ) =-A*COSH ( A* ( ETA +H ) ) *S l NH ( A* ( Y +H ) ) *COS 
1 ( A* ( X- X I ) )*A/ (C0SH(A*H)**2*(  l A* A- ANU* ANU ) *H  + ANU)  ) 

Q20Y Y ( Y , ETA, X, XI) =-A*A*CQSH(A*(Y  + H))*C3SH(A*( ETA  + H))* 
1C0S( A*(  X-X I ) )*A/(COSH(  A*H) **2*(  ( A* A- ANU* ANU ) *H  + ANU  ) ) 
Q1S(Y , ETA, XX, AMU) =-P 1 (Y,ETA,XX,AMU)/( AO+AMU*H*( AA+Bd* 

1 (AMU- A) +CC*( AMU-A)**2+DD*(  A MU- A )**3)  ) 

Q2X  S ( Y , ETA , X,X I , AMU )=-P2X(Y,ETA,X,Xl,AMU)/<  AO  + AMU*H* 

1 ( AA  + BB* ( AMU- A) +CC*( AMU-A)**2+DD*(AMJ-A)**3)  ) 

Q 2YS(  Y , E T A , X , X I , AMU ) =-P 2Y ( Y , E T A , X.XI , AMU) / ( AO  + AMU*H* 

1 (AA+BB* ( AMU-A) +CC*( AMU-A) **  2 + DD* ( AMU-A>**3) ) 

Q2YY  S ( Y , ETA, X.XI  , AMU) =-P2YY (Y ,ETA,X,X I , AMU J /( AO+AMU*H* 
1 ( AA  + BB* (AMU-A ) +GC*< AMU-A ) **2+0D* ( AMU-A  1**3)  ) 

FUN1 ( Y , ETA ,XX, AMU) =EXP< -AMU*H) *S INH< AMU*ETA)*SINH 
l ( AMU*Y ) *COS( AMU*XX) /(  A MU* COSH  ( AMU*H  ) ) 

FUN2(  Y , ETA, XX) = A . * 3 . IA159*C0SH( A*(Y+H) ) *CCSH( A* 
1(ETA+H))*CQS(A*XX)/(2.*A*H+SH2AH) 

FUNR ( X , Y , X l ,ETA , aNX, ANY ) = ( l X-X I ) * ANX  + ( Y + ET  A ) * ANY ) / 

1( (X-XI)**2+(Y+ETA)**2) 

FUN3X  ( Y »ETA,X,XI ,AMU)=EXP(  - A MU *H ) *S I NH  ( AMU*ETA)*S I NH 
l ( AMU* Y ) *S IN<  AMU* ( X-X I ) ) /COS  H ( AMU  *H ) 

FUN3Y ( Y ,ETA, X,XI ,AMU) =-EXP ( -AMU*H ) *S I NH ( AMU*E  T A ) *COSH 
1 ( AMU* Y ) *C OS ( AMU* ( X -X  I ) ) /C OS H(  A MU*H ) 

FUN3YY( Y ,ETA , X, XI , AMU) =-EX P < - AMU*H) *S INH( AMU* ETA) * 
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IS  INH< AMU*Y )*AMU*AMU*COS( AMU*(X-XI ) ) /COSH l AMU*H) 
FUN4(  Y,  ETA,  XX,  A, NIX,  ANY  .SIGN)  =4. *3.  14159*A*CCSH(A* 

1 ( ETA«-H)  ) *(  ANY*SINH(A*(Y+H  ) » *CGS  ( A*XX ) - A N X*C  OSH  ( A* 
2(Y+H))*SINlA*XX)*SIGN)/(2.*A*H-t-SH2AH> 


EVALUATION  OF  THE  FINITE  INTEGRAL  IN  G TO  DETERMINE 
THE  SIZE  OF  THE  SUBDIVISIONS 


XX=ABS( X-XI  ) 

IF(X.LT.XI)  S I GN  =-  1 . 0 
IF(X.GT.XI)  SIGN =1.0 
IF(X.EO.XI)  S I GN  = 0 . 0 
DO  50  N = 1 , 15 
DELMU=2*A/  ( o-*N  + 3 ) 

AMU=0. 0 
SUM=0.0 

F 0 = ( Q 1 ( Y, ETA, XX, AMU)-Q10(Y , ETA, XX) )/( AMU- A) 

LL  = 6*N+  3 
CO  40  NN=i,LL 

IFIABSI  AMU-A! .LT. .00001)  GO  TO  10 
IF (A3S(  AMU+DELMU-A) .LT. .00001 ) GO  TO  10 
IF  I A3S I AMU+DELMU/3.-A  ) .LT . .00001  ) GO  TO  10 
IF(A3S( AMU*2.*DELMU/3.-A) .LT..00001)  GO  TO  10 
F 1=(Q1  ( Y , ETA, XX,  AMU+DE  LMU/3 • ) ’■Q  10 ( Y ,ETA,  XX)  )/(AMU«- 
1DELMU/3 .-A) 

F 2=  ( Q1  ( Y,  ETA,  XX,  AMU+2 .*DELMU/3 . )-Q10[Y,ETA,XX))/ 

1(  AMU+2.*DELMU/3.-A) 

F3=(gi( Y, ETA, XX, AMU  + OE  LiMU  ) -Q 1 0 1 Y , ETA, XX) ) /( AMU+OELMU 
1- A3 

GO  TO  30 

10  F0=(Q1( Y, ETA, XX, AMU+DELMU ) -Q10 l Y , ETA, XX) ) /( AMU *0 EL MU 
1-A) 

GO  TO  40 

30  SUM=  (DELMU/3.  ) * ( FO +3 .*F 1+3. *F2 +F3 } + SUM 
F0  = F3 

40  AMU=AMU+DELMU 
TEST! N) = SUM 
IF(N-l)  50,50,45 
45  MN=N-1 

MM=6*MN+3 

I F ( AB  S ( (TEST(N)-TEST(N-l)  )/TEST(N)).LT.  .010)  GO  TO  <>C 
50  CONTINUE 
60  CONTINUE 
PV1F=2.*SUM 


EVALUATION  OF  THE  FINITE  INTEGRAL  IN  GN  USING 
2*A/MM  SUBDIVISION  SIZE 


DELMU=2*A/MM 
AMU=0. 0 
SUMX=0 . 0 
SUMY  =0.0 

F0= (Q2X ( Y,ETA,X, XI , AMU ) -Q20X ( Y , ETA,X,XI  ) ) /( AMU- A) 
Y0=(Q2Y(Y»ETA,X,XI ,AMU)-320Y(Y,ETA,X,XI ) ) / ( AMU- A) 

DO  90  N N = 1 , M M 

IFIABSl AMU-A) .LT. .00001)  GO  TO  70 

I F ( ABS ( AMU+DELMU-A) .LT. .00001)  GO  TO  70 

IF ( A6S< AMU+OELMU/3 .-A ) .Li.  . 00001 ) GO  TO  70 

IF( ABSl AMU+2.+3ELMU/3.-A) .LT. .00001)  GO  TO  70 

F 1 = ( Q 2 X ( Y , ETA, X, X I , AMU+DELMU/3 . ) -Q20X ( Y ,ETA,X,XI))/ 

I ( AMU+OELMU/3. -A) 

F2=(Q2X(Y  , ETA  , X , X I , AMU+2.  *DELMU/3  . ) -Q20X  ( Y , ETA  , X , X I ) )/ 
1 ( AMU+2. ODELMU/3 .-A) 

F3=IQ2X<Y,ETA,X,XI  , AMU  + DE  LMU  ) - Q2  CX  ( Y , E T A , X , X I ) )/ 

1 ( AMU+OELMU-A  ) 

Y1=(02Y(Y,ETA,X,XI , AMU+DELMU/3 . ) -C20Y (T,ETA,X,XI))/ 

1 ( AMU+DELMU/3. -A ) 

Y2=(Q2Y  ( Y,  ETA,X,X  l , AMU+2 .*0ELMU/3 . )-Q20Y  (Y,£TA,X,XI))/ 
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l(AMU*2.*DELMU/3.-A) 

Y3=(Q2Y(Y,ETA,X,Xl , AMU+DELMU) ^020Y (Y,ETA,X,XI  ))/ 
1 ( AMU*DELMU-A  ) 

GO  TO  75 
70  CONTINUE 

F 0= ( Q2X ( Y , ETA, X, XI , AMU +DELMU ) -Q2CX ( Y,ETA, X, XI ) } / 
1 ( AMU*DE  LMU-A  ) 

Y0= ( Q2Y ( Y , ETA, X,XI , AMU+DELMU ) -Q20Y (Y,ETA,X,XI  ) )/ 
I l AMU+DELMU-A) 

GO  TO  30 
75  CONTINUE 

SUMX=SUMX*  (DELMU/8. )*{F0+3.*F1+3.*F2+F3) 

SUM Y= SUMY*  (DELMU/8. )* ( YO +3 . * Y 1 *3 . *Y2 *Y3  ) 

F0  = F3 
Y 0 = Y 3 

80  AMU=AMU  *DELMU 
PV2FX=2.*SUMX 
P V2FY=2 .*SUMY 


EVALUATION  CF  THE  FINITE  INTEGRAL  IN  GY Y USING 
2 111  A/MM  SUBDIVISION  SIZE 


DELMU=2*A/MM 
AMU=0.0 
SUMYY  = 0.0 

YY0=(Q2YY(Y,ETA,X,XI  , AMU) -220Y Y ( Y ,ETA , X , X I ) J/IAMU-A) 

DO  350  NN  = 1 , MM 

I F ( ABS ( A MU- A ) .LT .0.00001)  GO  TO  320 
IF  (ABS( AMU*DELMU-A) .LT. 0.00001 ) GO  TO  320 
IF  (A3S(AMU*DELM'J/3. -A). LT. 0.00001)  GO  TO  320 
IF  (ABS(AMU*2.*DELMU/3.-A).LT.0.00001)  GO  TO  320 
YYL  = ( Q2YYI Y,ETA, X,XI , AMU  *DE  LMU/ 3 .)-Q23YYlY,ETA,X,XI))/ 
1 t AMU+DELMU/3.-A ) 

YY2=(Q2YY(Y,ETA»X,XI,AMU+2.*0ELMU/3.) -0 20YY ( Y , E T A , X , 

IX  I ) )/( AMU*2.*0ELMU/3.-A) 

YY3=(  Q2YY(  Y,  ETA,  X,XI  , AMU+-DELMU  ) -Q20YY  ( Y , ET  A , X , X I ) ) / 

1 ( AMU*DE  LMU-A  ) 

GO  TO  325 
320  CONTINUE 

YY0=(Q2YY(Y»ETA,X,XI, AMU*DE LMU ) -Q20YY (Y,ETA,X,XI))/ 

1 ( AMU+DE LMU-A ) 

GO  TO  350 
325  CONTINUE 

SUMYY  = SUMYY* ( DELMU/8. ) *(YY0*3. *YY1+3.*YY  2 + YY3) 

Y YO  = YY  3 

350  AMU  = AMU  + DE  LMU 
PV2FYY=2.*SUMYY 


EVALUATION  OF  THE  INFINITE  INTEGRAL  IN  G,  GN,  AND  GY Y 
SIMULTANEOUSLY 


AMU=2*A 
DELMUO=  DELMU 

F0=Q1  ( Y ,£TA,XX, AMU) /( AMU- A) 

F 0X=3  2X ( Y , E T A , X, XI ,AMU) /( AMU- A) 

F0Y=Q2Y ( Y ,ETA,X, XI , AMU) / 1 AMU- A) 

F0YY=Q2YY(Y,ETA,X,XI ,AMU) /( AMU- A) 

SUM=0 . 0 
SUMX=0.0 
SUM Y=0 . 0 
SUMYY  = 0 .0 
DO  100  NN= I ,200 
DO  95  N=1 , 20 

Fl=Wl  ( Y , ETA, XX, AMU*DELMU/3.  ) / ( A MU *DEL M J / 3 . - A ) 

F 1X  = Q2X(  Y ,ETA ,X , XI , AMU*DELMU/3 . ) / l AMU  *D  E LMU  / :> . - A ) 

F 1Y  = Q2Y ( v , ETA , X , X I , AMU *0E L MU / 3 . ) / ( AMU*OE LMU/3  . -A ) 
F1YY=Q2YY(Y,ETA,X,XI, AMU* DEL  MU/ 3 . ) / ( AMJ * DE L MU/ 3 . - A ) 
F2=Ql ( Y ,ETA, XX, A MU* 2. *DE LMU/3. ) / ( AMU* 2 . * DEL  MU/ 3 . - A ) 
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F2X=Q2X ( Y, ETA, X, XI , AMU+2 .*DELMU/3 . > / ( AMU* 2 . *DEL MU/ 3 . 

1 — A J 

F2Y=32Y(Y , ETA, X, XI , AMU  + 2 . *DE LMU/3 . ) / ( AMU  + 2 . *DELMU/ 3 . 

1- A) 

F2YY=Q2YY( Y,  ETA,  X,  XI  , AMU+2. *DE LMU/3. )/(  AMU  + 2.  *DELMU/ 3. 
I-A) 

F3  = Q1 { Y , ETA, XX , AMU*  DELMU )/ ( AMU+DELMJ-A) 

F3X=Q2X ( Y , ETA ,X ,XI , AMU+DELMU) / l AMU  + DEL  MU- A ) 

F3Y  = U2Y ( Y, ETA,X , X I , AMU+DELMU) / ( AMU+DELMU- A ) 

F3YY=Q2YY<  Y,ETA,X,XI, AMU  + DELMU) /( AMU  + DELMU- A) 

SUM  = (DELMU/3. )*(FO+3.*Fi+3.*F2+F3) +SUM 
S UMX=  (DELMU/8  .)*(F0X  + 3.*F1X  + 3.*F?X  + F3X)+SUMX 

SUMY=  ( DELMU/3. ) * ( FOY+3 . *F 1 Y +3 . *F2Y +F 3 Y ) + SUMY 

SUMYY  = SUMYY+( DEL  MU/8. )*(F0YY+3.*F1YY+3.*F2YY+F3YY) 

F 0 = F 3 
FOX=F3X 
FOY=F3Y 
F QYY  = F3 YY 

ASM=EXP(AMU*(ETA+Y) ) /AMU 
IF ( ASM. LT . .0001 ) GO  TO  86 

95  AMU= AMU+DELMU 
86  T E ST( NN ) = S UM 

TESTT(NN)=SUMX 
TESTTT ( NN ) =SUMY 
TESTYY ( NN  ) =SUMYY 
IF(NN-l)  100,100,97 

97  CONTINUE 

IF( ASM. LT. .00010 ) GO  TO  105 
IF ( ABSl SUM) .LT. .0001 ) GO  TO  98 

1F(A3S( ( TEST ( NN) -TEST (NN-1 ) ) / T t S T ( NN  ) ) . G T . . 00  1 ) 

1GG  TO  100 

98  IF (ABS( SUMX) .LT. .0001 ) GO  TO  99 

I F ( AB  S < ITESTTINN ) - TE S TT ( NN- 1 ) )/TESTT(NN)).GT..001) 

1 GO  TO  100 

99  IF  U6S(  SUMY)  .LT.  .0001  ) GO  TO  93 

I F(ABS(  ( TESTTT ( NN ) - TE ST T T ( NN- 1 ) ) / TE S TT T ( NN ) ) . GT . . 00 1 ) 
1G0  TO  100 

93  IF(ABS( SUMYYl.LT. 0.0001)  GO  TO  96 

IF  ( ABS ( ( TESTYY i NN ) -TESTYY I NN-1 ) ) /TESTYY (NN ) ) .GT. .001 ) 
1 GO  TO  100 

96  CONTINUE 
GO  TO  105 

100  CONTINUE 

102  WR I TE ( 6 , 1 0 3 ) 

103  FORMAT (3X3 4H INFINITE  INTEGRAL  DID  NOT  CONVERGE) 

105  CONTINUE 

PV1I=2.*SUM 
PV2IX=2 .-SUMX 
PV2IY=2.*SUMY 
PV2IYY=2.*SUMYY 
D ELMU= . 2/H 

IF  (J.EQ.25)  GO  TO  2A0 

AMU=0 . 0 

SUM=0. 0 

SUMX  = 0 . 0 

SUMY=0. 0 

SUMY Y = 0.0 

F0=0 .0 

F0X=FUN3X ( Y, ETA, X, XI, A MU) 

F0Y=FUN3Y( Y,ETA, X ,XI , AMU) 

F0YY  = FUN3YY(Y,ETA,X,XI  ,AMU) 

DO  200  NN  = 1 , 1 0 0 

IF ( AMU*H.GT.5. ) DELMU=.l 

DU  195  N= l , 20 

F1X  = FUN3X(Y,ETa,X,XI, AMU  + DELMU/ 3.  ) 

F 2X=FUN  3X ( Y , ETA, X, XI  , A MU+2. ^DELMU/3. ) 

F 3X  =FUN  3 X ( Y , ET  A , X , X I , AMU  + DELM'J  ) 

F 1 Y = FUN3 Y ( Y, ETA, X , XI , A MU  + DE LMU / 3 . ) 

F2Y=FUN3Y( Y,ETA,X,XI , AMU+2. *DE  LMU/3. ) 

F 3Y  =FUN  3 Y I Y,ETA,X,XI , AMU  + DELMU ) 

F 1YY  = FUN3YY( Y,ETA , X,X I , AMU+ DEL  MU/ 3 . ) 
F2YY=FUN3YYIY,ETA,X,XI , AMU +2.* DELMU/3. ) 
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F3YY=FUN3YY( Y,ETA,X,XI , AMU+DELMU) 

I F ( AMU  + DE  LMU .LT..001)  GO  TO  120 
F 1 = FUN1  (Y  »ET  A»XX,  AMU  + DE  LMU/ 3 . ) 

F 2 = FUN  1 (Y  fETAtXXjAMU  + Z.^DEL  MU/  3 . ) 

F3=FUN1 ( Y , ETA, XX , AMU  +DELMU ) 

GO  TO  130 
120  CONTINUE 

Fl=EXP(-( AMU+DELMU/3.  ) *H ) *C OS ( ( AMU+QELM J / 3 . )*XX)* 

1 ( AMU+DELMU/3 . ) *Y*E  T A/COSH ( ( AMU  + DELMU/3.  ) *K> 

F2=EXP(-(  AMU+2.*DELMU/3.  )*H)*CGS(  ( AMU + 2 . *DE  LMU/ 3 . )*XX) 
l*(AMU+2.*DELMU/3. ) *Y* ETA/ COSH ( ( AMU+2 . *DE LMU/3 . ) *H  ) 
F3=EXP(-(AMU+DELMU)#H)*C0S( ( AMU+DELMU ) *XX ) * ( AMU+DELMU) 
l *Y*ETA/COSH( ( AMU+DELMU ) *H ) 

130  CONTINUE 

SUM=  (DELMU/8. ) * ( FO +3 .*F 1 +3 . *F2 +F3 ) +S UM 
SUMX-  ( DEL  MU  / 8 . ) * ( FO X+3 . * F IX  + 3 .*F2X +F 3 X ) + SUMX 

SUMY=  (DELMU/8 . ) * ( FOY  + 3 . *F 1 Y + 3 . *F 2Y  + F 3 Y ) + SU MY 

SUMYY=SUMYY+(  DELMU/8.  )’MF0YY  + 3.*FIYY  + 3.*F2YY  + F3YY) 
F0=F3 
F 0X=F3X 
F0Y=F  3Y 
FQYY=F3 YY 

ASM=EXP( AMU*(ETA+Y)  ) 

IF(ASM.LT.. 0001 ) GO  TO  196 

195  AMU= AMU  +DELMU 

196  TE  ST ( NN ) = SUM 
TESTT(NN)=SUMX 
TESTTT( NN) =SUMY 
TESTYY ( NN) =SUMYY 
IF(NN-l)  200,200,199 

199  CONTINUE 

IF (ASM. LT. .00010)  GO  TO  205 
IF(ABS( SUM) .LT. .0001  ) GO  TG  206 

IF(ABS( ( TEST(NN) -TEST (NN-1) )/TEST(NN) ) .GT. .0010) 

1 GO  TO  200 

206  I F ( ABS ( SUMX) .LT. .0001  ) GO  TO  207 

IF  (A8SI  ( TESTT(NN  ) — TES-TT  ( NN-  l))/TESTT(NN)).GT..0010) 

1 GO  TO  200 

207  IF ( ABS ( SUMY) .LT . .0001 ) GO  TO  204 

IF (ABS ( (TESTTT(NN)-TESTTT(NN-l) ) / TEST T T ( NN ) ) . GT . . 00 1 0) 
1 GO  TO  200 

204  I F < AS  S ( SUMYYJ.LT. 0.0001)  GO  TO  208 

IF  (A8S( ( TESTYYl  NN)-TESTYY(  NN-1 ) ) / ( TEST YY ( NN) ) ) . 

1 GT. .001 ) GO  TO  200 

208  CONTINUE 
GO  TO  205 

200  CONTINUE 
WRITE(6 ,202) 

202  FORMAT ( 3X14HN0  CONVERGENCE) 

205  CONTINUE 
GINF=-2*SUM 
GXINF=2 .*SUMX 
GY  I NF  = 2 . *SUMY 
GYYINF=2.*SUMYY 
GNSING=  .5 

IF  (I.E0.25)  GO  TO  218 
I F( I . EQ . J ) GO  TO  220 
A I J J-  I - J 

THETA1= ABS< A I JJ) *D  EL  THE 

A I N J J = I +NPTS-J 

T HET A2  = ABS ( A IN J J ) *DEL THE 

AJNII=I-J-NPTS 

THETA3  = ABS( AJNI I )* DEL  THE 

THETA=THET  A1 

IF(THETA2.LT. THETA)  THETA=THETA2 
IF(THETA3.LT. THETA)  THET  A = THET  A 3 
IF ( THETA. GT.. 15)  GO  TO  218 

GS I NG=D EL THE *ALOG( 2 . ) +2 . * ( -DEL  THE / 2. * ( DEL TH E/ 4 . + 

1 THETA/ 2 . )*ALOG<  T HE T A/ 2 . + OEL T HE/ 4 . ) - ( THE  T A/ 2 .-DEL  T HE  / 4. 
2) *ALOG( THETA/2. -DEL THE/ 4. )-(0ELTHE/4.+THETA/2. ) **3/ IS. 
3+ (THET A/ 2. -DEL  THE/4. )**3/18.-( THETA/2. +DELTHE/ '+ . )**5/ 
4( 180. #5. ) + ( THETA/2. -OELTHE/4. ) * * 5 / ( 180. *5. ) -( THET  A / 2 . + 
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5DELTHE/4. ) **7/(  2 83 5.  *7.  ) + ( THE T A / 2. -DELTHE / 4 . ) ** 7/ 
6(2835. *7.  ) ) 

GXSING= ( X-XI )/( ( X-XI )**2M Y-ETA) **2) 

GYSING=( Y-ETA) /(  (X-XI )**2+(  Y-ETA )**2) 

GO  TO  230 

218  GSING=DELTHE*ALOG(SQRT(  ( X-X  I ) **2  + l Y-ETA  ) *<=2  ) ) 

GX  S i NG  = (X-XI >/ ( ( X-XI )**2+( Y-ETA )**2) 

GY SING= (Y-ETA) /( ( X-XI ) **2 + ( Y-E T A ) **2 ) 

GO  TO  230 
220  CONTINUE 

GSING=DELTHE*AL0G(0ELTHE/2. )-OELTHE- (DELTHE**3)/ 
l(18.*16.)-(DELTHE**5)/(  1 80 . *5 . * 2 56 . ) - ( 3 E LTH E** 7 ) / 
2(2835. *7. *256. *16. ) 

GXSING=GNSING*ANX 
GYSING=GNSING*ANY 
230  CONTINUE 

G =GSI  NG/DELTHE-ALOG( SQRT( ( X-X I ) **2  + ( Y * ET A )**2 ) ) + P V IF 
UPVIH-GINF+CMPlX  (0.  » -1.  ) *FUN2(  Y , ETA, XX  ) 

G X=GX  S I NG-  FUNP  (X,  Y,  XI,  ETA,  1.0,0. 0) +PV2FX*  PV 2 1 X *GX I NF  + 

1 F UN4(  Y , ETA ,XX, 1 .0 , 0.0 , S IGN) *CMPLX (0.0, -l .0) 

GY  = GY  S I NG-FUNR( X » Y ,X  I , ETA ,0 .0* 1 .0 ) +PV2FY  +PV2IY+GY I NF  + 
IF UN4(Y, ETA, XX,  0.0 ,1.0, SIGN) *C  M PLX (0.0»  — 1.0) 

IF  (I.LE.NPTS)  GO  TO  235 

FUPYY=1./SQRT(((X-XI)**2+(Y-ETA)**2)**3)-((Y-ETA)**2+ 
1( Y-ETA) )/SQRT( ( (X-XI )**2M Y-ETA) **2)**5) 

FUNYY=I ./SQRT( ( ( X-X I ) * *2  + ( Y + E T A ) **2 ) ) - ( < Y+ ET A ) **2  + 

1 ( Y + ETA)  ) / S QRT  ( ( ( X- X I ) **  2 *■  ( Y *E  T A ) **  2 ) **  5 ) 
GYY=FURYY*-FUNYY  + PV2FYY*PV2  I YY<- GY Y INF  + G 'IP LX  ( 0 . 0 , - 1 . 0 ) * 

1 FUN2 ( Y ,ETA,XX)*A*A 
235  CONTINUE 

IF  (I.EQ.25)  GO  TO  250 

GN=GNS I NG-f JNR { X » Y ,X  I , ET A , A NX , ANY )♦ ( P V 2 F X + P V2 I X + G X I NF ) 
l*ANX+(  PV2FY+PV2IY+GYINF)*ANY*FUN4(Y,ETA,XX  , ANX,  ANY, 
2SIGN)*CMPLX(0.0,-1.0) 

GO  TC  250 
240  CONTINUE 

G = - ( P V 1 F + P V 1 1 ) +CHPLX(0. ,-l. )*FUN2( Y, ETA, XX  ) 

GN  = - ( PV2FX+PV2IX) *ANX- ( PV2F  Y *P V 2 I Y ) *ANY+FUN4( Y, ETA, XX, 
1 ANX, A NY  »S!GN)*CMPLX(0.0,-1.0) 

250  CONTINUE 
RETURN 
END 


THIS  SUBROUTINE  CALCULATES  G AND  ITS  DERIVATIVES, 

GX,  GY , GN , AND  3YY  , BY  USE  OF  THE  SERIES  FORM  FOP.  THE 
CASE  (X(I)  - XIJ))  GREATER  THAN  SHIN 


SUBROUTINE  GREENS  ( A , ANU, SH2 AH , SHAH ,CHAH , CQSA MU , S I N AMU 

1 , AMU, COE FG , SHY  I ,CHYI,SHYJ,CHYJ,I,J,XV,YV,XC,YC,ANXI, 

2 ANY  I , AN  X J , ANY  J , G I J , GN I J , GN J I , *X I J ,GX J I , G Y I J , GY J I , G YY ) 

COMPLEX  GI J,GNI J ,GNJI ,GXI J, GXJ I , GY  I J,GY J I , GYY 
DIMENSION  COS A MU ( 200, 25) , SINAMU( 200,25) , AMU (2 00) , 
1C0EFG(  200) 

DIMENSION  TESTK40)  ,TEST2(40)  ,TEST3(40)  ,TEST3G(40)  , 

1 TE  ST4( HO) 

COMMON/ CONST /H, J , DE LTHE , SMI N ,NPTS  ,NSPTS 

CCMMON/C  P I / P I 

SUM  1 = 0 . 0 

SUM2=0 . 0 

SUM3=0. 0 

SUM30=0.  J 

SUM4=0. 0 

DO  20  N = 1 , 40 

DC  15  K K= 1 ,5 

K=(N-1 ) *5+RK 

SNI =S INAMU ( K , I ) 

SN J=S I NAMU ( X , J ) 

CSI=COS AMU ( K , I ) 

CSJ=COSAMU( K , J ) 
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AK=AMU ( K) 

CCK  = COE  FG I K ) 

VAL=AK*ABS(XV-XC  ) 

IF  (VAL.GE. 75.0)  GO  TO  10 
EXI J=EXP( -VAL  ) 

GO  TO  11 

10  EX l J = 0 . 0 

11  CONTINUE 

SUM  1 = SU  M l+C  OK  *C  S I *C  S J *E  X I J 
SUM2=SUM2+-COK*CSI*CSJ*EXIJ*AK 
SUM3  = SUM3*COK.*AK*SNI  *CSJ*EX  I J 
SUM3G=SUM3Q*CQK*AK*SNJ*CSI*EXIJ 
SUM4=SUM4+-AK.*AK*C0K*CSI*CSJ5!‘EXIJ 

15  CONTINUE 
TESTl ( N ) =SUM1 
TEST2I N ) =SUM2 
TEST3I N) = SUM3 
TEST30(N)=SUM30 
T ES  T4 1 N ) = S UM4 
IF(N-l)  20,20,16 

16  IFIABSI  SUM1) .LE. 0.000001)  GO  TO  17 

IF(  ABSI  ( TEST  1(NI  -TEST  l(N-l)  ) / TEST  1 1 N ) ) . GT . . 00  10 ) GO 
1 TO  20 

17  I FIABSI  SUM30) .LE .0.000001 ) GO  TO  13 

IFIABSI ITE$T3QIN)-TEST3QIN-1) J/TEST3QIN) ) .GT. .0010) 

1 G 0 TO  2 0 

18  IFIABSI SUM3) .LE. 0.000001)  GO  TO  19 

IFIABSI  I TEST3I  N) -TEST31N-1 ) ) / TE ST3 I N) ) . GT . . 00 10 ) GO 
1T0  20 

19  IFIABSI SUM2) .LE. 0.000001)  GO  TO  21 

IFIABSI < TEST2I NJ-TEST2 IN- 1) 1/TEST2IMI J . GT.. 0010)  GO 
1T0  20 

21  IFIABSI SUM A ) .LE. 0.000001)  GO  TO  30 

IF! ABS I I TEST 4 1 NJ-TEST4 I N-l ) J/TEST4IN) ) .GT. .0010)  GO 
1T0  20 
GO  TO  30 

20  CONTINUE 
WRITEI6.25)  I.J 

25  FORMAT  I 10X23 HG REE NS  DID  NOT  CONV ERGE , 2X 2H I = I 2 , 

12X2HJ= 12) 

30  CONTINUE 

IF  IXV.LT.XC)  S I GN  = - l . 0 
IF  (XV.GT.XC)  S I GN  = l . 0 
IF  (XV.EQ.XC)  S I GN=0 . 0 

TERM4=4.*PI*CHYI *CHYJ*SIN{ A*A6S I XV-XC ) ) / I 2.  *A*H<-SH2AH) 
TERM5=4.*PI*CHYI  *CHYJ*CCS<  A* ABS  I XV-XC  ) ) / I 2.  *A*H-t-SH2AH) 
TERM6=A*TERM5*SIGN 
TE  RM7=A  *TE  RM4#  SIGN 

TER  MB  =4.*PI*A*SrlYI*CHYJ*C0S(A*ABSI  XV-XC  ) ) /(  2.*A*m- 
1SH2AH) 

TERM9=4.*PI*A*SHY I*CHYJ*SINI A* ABSI XV-XC  ) ) /( 2.*A*H  + 
1SH2AH) 

TERM80=4.*PI*A*SHYJ*CHYI*C0S(  A* ABSI XV-XC  ) ) / <2.*A*H* 
ISH2AH) 

TERM9Q  = 4.*PI$A*SHYJ<:CHYI*SIN(A*ABSIXV-XC  ) )/(2.*A*H+ 
1SH2AH) 

TERM10=A*A*TERM* 

TERM11=A*A*TERM5 

IF  IJ.EQ.25)  GO  TO  70 

G I J=CMP  LXt SUM1  + TERM4,-TERM5 ) 

GXI J=  CMPLXI -SIGN*SUM2*TERM6 ,TERM7) 

GY  I J=CMPLX<-SUM3+TERM9,-TERM8) 

IF  I I.EQ.25)  GO  TO  40 
GXJI=CMPLX( SIGN*SUM2-TERM6,-TERM7) 

GYJI=CMPLX( -SUM30+TERM90,-TERM8Q) 

GNI J = CMPLX l-ANX I*S IGN*SUM2*ANX I *TERM6-ANY I * SUM3*ANYl* 
1TERM9,ANXI*TERM7-ANYI*TERM8 ) 

GNJI  = CMPLXI ANXJ  *SI GN*SUM2-ANX J*TERMo-ANY J*SUM30+ANY J* 
1TERM90, -ANXJ*TERM7-ANYJ*TERM80) 

GO  TO  60 
40  CONTINUE 

G YY=CMP LX  I -SJMA+TERM10 ,-TERMl 1 ) 
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60  CONTINUE 
GO  TO  SO 

70  GI  J=CMPLX(-SUM1-TERM4,-TERM5) 

GNI J=CMPLX ( ANXI*S I GN* SUM2- A NX  I * T ERM6  + AN Y I *SUM3-ANY I * 
1TERM9,ANXI*TERM7-ANYI*TERM8  ) 

80  CONTINUE 
RETURN 
END 


THIS  SUBROUTINE  INVERTS  COMPLEX  MATRICES  TO  SOLVE 
THE  MATRIX  EQUATION  ALPHA (I , J ) *F ( I , 1 ) = HH (1,1)  AND 
ALPHA (IfJ)*Fl(IflJ  = PK ( I , 1 ) 


SUBROUTINE  COMA T ( N ,M  , A , B , D , I ) 

INTEGER  C fH, R, Q, Z 

COMPLEX  A,3,D,TT,P 

DIMENSION  A<N,N) ,B(N,M) ,C( 100,3) 

D = ( 1 .0, 0.0) 

DO  20  J = 1 , N 

20  C<  J ,3)  = 0 

DO  21  K = 1 , N 
TT  = (0.0,0.0) 

T = 0.0 
DC  4 J = 1,N 

IF  ( C ( J , 3 ) .EQ.  1 ) GO  TO  4 

DO  5 H = 1 ,N 

IF  ( C ( H , 3 ) - 1)  15,5,12 

15  IF  (T  .GE.  C ABS ( A( J , H ) ) ) GO  TO  5 

R = J 

Q = H 

T = CABS( A ( J , H)  ) 

5 CONTINUE 
4 CONTINUE 

C (Q,3)  = C<  Q, 3 ) * 1 
C ( K,  1 ) = R 
C(K,2)  = Q 

IF  (R  . EQ.  Q)  GO  TO  11 
D = -D 

DO  8 L = 1 , N 
TT  = A ( R , L ) 

A ( R ,L ) = A ( Q , L ) 

8 A ( Q , L)  = TT 

I F ( M .LE.  0)  30  TO  11 
DO  2 L = 1 »M 
TT  = B ( R , L ) 

B(R,L)  = 3 ( Q , L ) 

2 B(Q,L)  = TT 
11  P = A(Q,Q) 

A ( Q ,Q)  = ( 1.0, 0.0) 

DO  13  L = 1 , N 
13  A ( Q , L ) = A ( 0 , L ) / P 

IF  (M  .LE.  0)  GO  TO  1 
DO  3 L = 1,M 

3 BIQ,L)  = B ( Q , L ) / P 
1 DO  21  Z = 1 , N 

IF  <Z  . EQ.  Q)  GO  TO  21 
TT  = A ( Z , Q ) 

A(Z,Q)  = (0. 0,0.0) 

DO  16  L = 1 , N 

16  A ( Z , L ) = A ( Z , L ) - A ( Q , L ) *T  T 
IF  (M  . LE.  0)  GO  TO  21 

DO  17  L = 1 , M 

17  B(Z,L)  = 3 ( Z , L ) - B ( Q , L ) * T T 

21  CONTINUE 

DO  19  II  = 1 , N 
L = N > 1 - 1 1 

I F ( C ( L , 1)  .EQ.  C(  L , 2 ) ) GO  TO  19 
R = C ( L , 1 ) 

Q = C ( L , 2 ) 
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DO  7 K = l.N 
TT  = A ( K |R  ) 

A ( K , R)  = A(K,Q) 
7 A l K,Q)  = TT 
19  CONTINUE 

DO  18  K = 1 f N 
IF  (C(Kr3>  .NE. 
18  CONTINUE 
I = 1 
50  RETURN 
12  I = 2 
GO  TO  50 
END 


1)  GO  TO  12 


//GO.SV SIN 


DD  * 
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